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ABSTRACT 



A general combustion and heat transfer model for porous 
media subject to Darcy flow is formulated. The transient 
one dimensional model treats the combustion process in two 
phases. During the initial phase, combustion occurs within 
the porous medium. The second phase occurs when the exo- 
thermic reaction moves to the air inlet surface of the medium 
resulting in surface recession. The temperature dependency 
of the system parameters and thermophysical properties is 
taken into account. An analysis of combustion in a carbon 
porous medium is presented, as well as an assessment of the 
accuracy of the model. 
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I. 



INTRODUCTION 



Studies were conducted at the Naval Weapons Center, China 
Lake, California to assess the survivability of the F-18 
aircraft when exposed to open pool fires. A portion of the 
F-18 aircraft is constructed of graphite-epoxy materials. 

When this composite material is exposed to a severe thermal 
environment, the epoxy burns off at a low temperature (about 
400-500 degrees Fahrenheit) leaving behind a porous graphite 
mat. Experiments were conducted by J. Fontenot [1] to evaluate 
the combustion behavior of the porous graphite or carbon mat. 
Due to the complex behavior of the combustion observed for the 
carbon fibers, a mathematical model to simulate this behavior 
was developed by Vatikiotis [2] . The model was formulated to 
consider a porous medium composed of cylindrically shaped 
fibers, typical of fiber reinforced composite materials. The 
details of the initial effort were presented by Vatikiotis. 
Significant improvements have been made to the initial combus- 
tion model, and the model has been extended to consider 
spherically shaped particles, typical of thermal flow reactors. 
The initial work of Vatikiotis [2] , where it applies to the 
present investigation is included in the discussion for the 
sake of completeness. Other analytical investigations taking 
a similar approach, but concerned with other aspects of com- 
bustion or thermally active porous media, have been performed 
by Kordylewski [3] , Sahota and Pagni [4] , and Mehta, Sams and 
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Luss [5]. Kordylewski performed a numerical analysis of homo- 
geneous combustion in a two dimensional reactor. The objec- 
tive of Kordylewski ' s investigation was to show that the flow 
within the reactor influenced the conditions necessary for 
ignition. Sahota and Pagni considered the transient behavior 
of heat transfer in a porous medium when exposed to a fire. 
Their analysis was part of an overall investigation for deter- 
mining the internal stresses in a structure caused by a fire. 
Combustion of the porous medium was not a consideration. 

Mehta, Sams, and Luss investigated the transient temperature 
response of a packed-bed reactor when the entering fluid 
temperature was decreased. Their model assumed that the 
temperatures of the porous solid and the fluid were the same. 

Slattery [6], Yaron [7], and Kassoy [8] present three 
approaches for mathematically modelling a porous media. 
Slattery's and Yaron ' s approaches were similar in using inte- 
gral methods in their development. Slattery considered only 
mass transfer whereas Yaron treated transport processes in 
general. Kassoy' s approach was based on performing Taylor 
series expansions of energy transport processes with respect 
to a differential volume of a porous medium. An energy balance 
yielded two differential heat transfer equations, one for the 
porous solid and one for the fluid. A fundamental requirement 
of this approach is that the particle diameters remain small 
with respect to the thickness of the porous medium. 

Although the above investigations were concerned with 
specific aspects of thermally active porous media, the present 
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combustion model was developed to treat a larger class of 
problems (i.e., combined mass transfer, heat transfer, and 
combustion). In contrast to other investigations, the present 
model takes into account the temperature dependency of all 
the thermophysical properties of the air and the geometric 
parameters of the porous medium. The physical parameters such 
as porosity, permeability, and thickness affect the magnitude 
and behavior of pressure-driven air flow within the porous 
medium. As discussed by Kordylewski [3] , the air flow influ- 
ences the conditions for ignition. Therefore, sensitivity 
analyses were performed using the combustion model showing 
the effects of these parameters on combustion. In addition, 
a functional relation between the internal pore velocity and 
an initial uniform temperature distribution needed to sustain 
combustion is presented. The initial work of Vatikiotis [2] 
showed that the initial conditions are important in deter- 
mining whether or not sustained combustion will occur. An 
analysis is presented showing the dependence of combustion 
on the shape of the initial temperature. In addition, the 
effects of the boundary conditions are examined using the 
combustion model. Specifically, there is a change in the 
behavior of combustion when insulated boundaries on the porous 
solid are changed to permit heat transfer by radiation. The 
Semenov model is used to explain the behavior observed during 
the analyses discussed above. Although the Semenov model 
was originally proposed for homogeneous combustion, the 
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results 

geneous 



show that the Semenov model can be applied to hetero- 
combustion in porous media. 
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II. DESCRIPTION OF THE PROBLEM 



The problem under investigation is that of a porous slab, 
with a pressure differential across its thickness, in a severe 
thermal environment. The effect of the pressure differential 
is to induce air flow through the porous medium. This internal 
air flow produces two opposing effects, (1) internal convection 
heat transfer, and (2) a supply of oxygen for heat generation 
by combustion. Whether the porous medium moves towards sus- 
tained combustion or extinguishment depends on the interaction 
of these effects. A mathematical model was formulated to 
provide an understanding of this interaction and its effect 
on thermal behavior . 

The mathematical model was formulated as follows. Energy 
balances on the porous solid and the convected air provide 
heat transfer equations for each. The heat transfer mechan- 
isms included in the model are (1) conduction, (2) convection, 
and (3) radiation. In addition, nonvolatile combustion is 
included in the porous solid heat transfer equation as a heat 
generation term of Arrhenius type . 

The conservation of species law was applied to the oxygen 
molecule concentration. The resulting mass transfer equation 
includes the transfer mechanisms of (1) molecular diffusion, 
and (2) species transport by convection. A term accounting 
for the consumption of oxygen due to combustion is included 
as well. Darcy's law was taken as the constitutive equation 
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for flow through the porous medium. The pressure gradient, 
needed to calculate pore velocity, is obtained by solving a 
combined Darcy's law and continuity equation. 

As the behavior of the system occurs over a large tempera- 
ture range, all thermophysical properties which depend upon 
temperature , such as conductivity, viscosity, and air density, 
are treated as temperature dependent properties. Moreover, 
as combustion consumes the porous solid, changes in porosity 
and permeability are accounted for in the transient analysis. 

In order to provide a general model, a number of boundary 
conditions are included in the formulation of the model . A 
detailed discussion of boundary conditions is presented in 
Section III.H. The final system of equations to be solved are 
four transient, nonlinear, coupled partial differential equa- 
tions. The three heat and mass transfer equations were solved 
by a Galerkin formulation of the finite element method. The 
remaining combined Darcy's law-continuity equation was solved 
by a shooting method. The details of the solution procedure 
are presented in Appendix C. 
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III. THEORY AND BACKGROUND 



A. DESCRIPTION OF THE POROUS MEDIUM 

For this investigation, a porous medium is defined as a 
solid containing interconnected pores (i.e., path for airflow). 
There are two classes of porous media, consolidated and uncon- 
solidated. Consolidated porous media are those where the 
solid constituent or matrix remains rigid. Consolidated por- 
ous media may be comprised of either individual particles 
attached to one another (e.g., sintered metal), or a solid 
where a portion has been removed (e.g., charcoal). Unconsoli- 
dated porous media are comprised of discrete particles (e.g., 
granular beds) . The properties of an unconsolidated porous 
medium can change if the porous medium is agitated or settling 
takes place. Both consolidated and unconsolidated porous 
media can have spatially varying properties . 

The physical properties which characterize porous media 
(i.e., allow comparison of porous media without regard to class 

or form) are (1) porosity, (2) specific internal area, (3) pore 

size or diameter, and (4) tortuosity. These properties are 
common to both consolidated and unconsolidated porous media. 
Porosity, p, is defined as the ratio of void volume to total 

volume. The specific internal area, z, is the ratio of in- 

ternal surface area to bulk volume. The dimension is the 
reciprocal of length. Pore diameter is used to describe the 
pore system of a porous medium. The actual shape and size of 
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a pore may be irregular and complex, and would be difficult 
to describe geometrically. In practice, there are many con- 
ventions used to define the pore diameter. For this inves- 
tigation, a hydraulic diameter theory is used and is discussed 
later in this section. The tortuosity of a porous medium is 
defined as the ratio of the length of the flow path of a fluid 
particle to the straight line distance. Originally, tortu- 
osity was considered a kinematic property, but was subsequently 
adopted for characterizing porous media. Tortuosity is also 
discussed later in this section. Scheidegger [9] discusses 
various methods used to measure properties of porous media. 
Though most methods discussed are based on experimental tech- 
niques, the above properties for this investigation are calcu- 
lated based on an idealization of the geometry of a porous 
medium. 

The porous medium was modelled as shown in Figure III.l. 

In Figure III.l, D is the distance between the fiber or parti- 
cle centers, and d is the particle diameter. The geometry is 
idealized since the actual porous medium may be irregular in 
particle diameter and distribution. For cylindrical fibers, 
the porosity, which was defined as the ratio of void per unit 
volume is , 



p = 1 - j(d/D) 2 

and for spherical particles. 



(III.l) 
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p 



1 - £(d/D) 3 



(III . 2) 



The pore diameter is obtained by an expression proposed by 
Carman [10 ] , 



where z is the specific internal area or the particle surface 
area per unit volume of porous medium. Expression III. 3 is 
analogous to the more familiar form of mean hydraulic diam- 
eter, 4V/A, where V is the void volume and A is the wetted 
surface area. The specific internal area, z, based on the 
idealized geometry of Figure III.l is calculated by. 



for spherical particles. Expressions III. 4 and III. 5 assume 
one half the total internal surface area is effective for 
convection heat transfer. This fractional amount of total 
area was an estimate based on Fontenot's experimental results 
and does not apply to porous media in general. A more general 
approach, as presented by Scheidegger [9] , is the Kozeny 



<5 = 4p/z 



(III. 3) 



z 



JTT d/D 2 



(III. 4) 



for cylindrical fibers, and by 



z 




(III. 5) 
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equation given by. 



z = (sp 3 /m) 1/2 (III. 6) 

In equation III. 6, m is the specific permeability or hydraulic 
conductivity (discussed in Section III .B) , and s is the Kozeny 
constant. The Kozeny constant is usually taken as .2 for 
specific surface calculations. Advantages to using the Kozeny 
equation are (1) the calculated values of specific internal 
area are in fair agreement with the experimental values, and 
(2) the calculation is independent of particle shape. A 
disadvantage is that the Kozeny equation fails for highly 
porous fibrous media. The tortuosity, t , is the ratio of the 
length of the flow path for a fluid particle to the straight 
line distance. For the geometric configuration of Figure III.l, 
the tortuosity depends on the ratio d/D. Carman [10] pre- 
sents a table of measured tortuosity factors of various materials 
and geometries, and points out the differences between the 
analytical determinations of tortuosity. Here, we adopt his 
recommended value of 1.4. Since particle size decreases as 
the carbon is consumed, all geometric properties which depend 
on fiber or particle diameter are functions of time and posi- 
tion. The model assumes that the carbon matrix remains rigid 
as the particle diameter decreases, and thus, porosity in- 
creases with combustion. In his discussion on the packing 
theory of spheres, Scheidegger [9] states the highest reported 
porosity for spherical particles in a stable configuration is 
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.875. Carman [10] reports of investigations of fibrous porous 
media with porosities as high as .99. The change in particle 
diameter is accounted for, and will be discussed in the sec- 
tion on the Arrhenius expression (Section III.D) for carbon 
combustion. 

B. DARCY'S LAW AND PORE VELOCITY 

Reynolds number for porous media is defined by 

Re = (p ud)/y (III. 7) 

3 . 

where u is the local pore velocity, p is the mass density of 
air, and y is the dynamic viscosity. Depending on the magni- 
tude of the Reynolds number, the motion of the fluid may be 
dominated by molecular, viscous, or inertial effects. Most 
investigations of fluid flow in porous media have been in the 
range of Reynolds number where the flow is dominated by vis- 
cous and inertial effects. The Navier-Stokes equations are 
applicable in describing the fluid motion in this range of 
Reynolds number. However, because of the convoluted geometry 
and the necessary boundary conditions (i.e., u = 0 at a solid- 
fluid interface) , solution of the Navier-Stokes equations are 
difficult for a porous medium. Extensive experimental work, 
as discussed by Scheidegger [9] has shown that fluid flow in 
porous media is governed by Darcy's law for the range of 
Reynolds number where viscous effects dominate. The upper 
limit (high velocity end) Reynolds number reported in the 
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experimental investigations varied from .1 to 75. As Reynolds 
number increases, inertial effects increase. However, Boffa 
[11] has shown that for a fixed Reynolds number, inertial 
effects diminish with increasing air temperature. For the 
problems covered in this investigation, the Reynolds number 
did not exceed a value of 5 . 

Neglecting body forces, Darcy* s law for one dimensional 
flow is. 



m.dP. 
y dx 



(III. 8) 



where Q is the filter velocity or the volumetric flow rate per 
unit of cross-sectional area, and dP/dx is the pressure gradi- 
ent. The specific permeability or the hydraulic conductivity 
of the porous medium is defined by. 



m = 96 p (6 /t) 



(III. 9) 



Expression III. 9 is based on a capillaric-serial model given 
by Scheidegger [9], It should be noted that, in a physical 
sense, permeability and porosity are not related. Porosity 
is a quantifiable property of the porous medium. Whereas, 
permeability is a constitutive property (i.e., a property 
specified by a constitutive equation, Darcy's law). As in 
the case of Fourier's law of heat transfer, Darcy's law is 
not derived from first principles, but has been obtained 
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through exhaustive experimental analyses. The Dupuit- 
Forcheimer assumption, presented by Carman [10] , relates 
the local pore velocity to the filter velocity by. 



Q - pu (III. 10) 

The hypothesis of the Dupuit-Forcheimer assumption is that the 
local pore velocity is greater than the filter velocity. 

Noting that the actual velocity in a single pore is not a 
constant, but rather a function of the location within the 
pore, the Dupuit-Forcheimer assumption defines an "average" 
velocity within the pore. 

The continuity equation for one dimensional fluid flow in 
porous media with nonconstant porosity distribution is, 

D(pp a ) 3u 

Dt~~ + pp a ^ = 0 (III. ID 



Substituting in the Dupuit-Forcheimer assumption and Darcy's 
Law, the continuity equation becomes. 




,_1_ ^a_ + 1 3m _ 1 3y,dP 
3x m 3x y 3x^ dx 

a 



u 

p a m 



3 (pp_) 

a 

3t 



0 (III. 12) 



Equation III. 12 along with the boundary conditions can be 
integrated to provide the pressure and the pressure gradient 
distributions through the porous medium. The boundary condi- 
tions are P(0) = P , the ambient pressure, and P (L) = P . The 
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pressure differential across the porous medium is AP = P^ - P^ 
Knowing the pressure gradient distribution, the pore velocity 
distribution is obtained from Darcy's law and the Dupuit- 
Forcheimer assumption. Derivations of equation III. 12 and 
the continuity equation are presented in Appendix A and B, 
respectively . 

C. SEMENOV MODEL OF COMBUSTION 

In this work we adopt the combustion model of N. N. Semenov 

as described in the texts of Frank- Kamenet ski i [12] and Vulis 

[13] . A brief discussion of those features of the model which 

relate to the present investigation will be given. Fundamental 

to the model is the relation of reaction rate to temperature, 

and the interaction of heat generation and heat transfer. 

The reaction rate, R , is represented by the Arrhenius law for 

c 

a simple n-th order reaction, 

R = A <]> n exp {——) (III. 13) 

c R T 

u 

where A ^ is the characteristic time of the chemical reaction, 

E is the activation energy, R is the universal gas constant, 

/v 

T is the absolute temperature, and <j> is the oxygen concentra- 
tion. A simple reaction is one in which the rate depends on 
the concentrations of the reactants and not on the products. 

The heat generated by the exothermal reaction is obtained by 
multiplying R c by the heat of combustion. 
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In Figure III. 2, heat generation R^ is plotted versus 
temperature. The resulting curve is referred to as the S- 
curve due to its sinuous shape. The S-curves are distinguished 
by two regions. In region I, the reaction rate and tempera- 
ture are low and as a result, there is an excess supply of 
oxygen. In this region the reaction is controlled by the 
temperature, and is called the kinetic regime of combustion. 

In the kinetic regime, the reaction rate increases exponen- 
tially with increasing temperature. In region II, the higher 
reaction rates are only slightly affected by the higher tem- 
peratures, and the reaction is limited by the availability of 
oxygen. Region II is called the diffusion regime of combus- 
tion. It should be noted that S-curves are not obtained by 
plotting Rg versus temperature for a constant oxygen concen- 
tration. The flattening of the S-curve at the higher tempera- 
tures is due to decreasing oxygen concentration at those 
temperatures. The S-curve, which represents the behavior of 
an actual system, demonstrates the strong coupling between 
temperature and oxygen concentration. 

In addition to the S-curves, Figure III. 2 shows two con- 
vection heat transfer curves, and / at air flow tem- 

peratures, T^ and T respectively. The equation for the 
heat transfer lines is of the form, 

q z = h(T s - T a ) (III. 14) 
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FIGURE III. 2 Semenov model of combustion. 



where h is a convection heat transfer coefficient, and T and 

o 

T are the solid surface and air temperatures, respectively. 

cl 

The intersection of the heat generation and heat transfer 
curves, = q^, defines the stable quasi-stationary point A. 
For a system in thermal equilibrium, the solid and air tempera- 
tures will remain at T A and , respectively. Perturbing the 
system from one equilibrium position to another causes the 
quasi-stationary point to adjust accordingly (e.g., raising 
the air flow temperature from to will cause the quasi- 
stationary point to move from T a to T t ) . At point I, the R 
and q^ curves are tangent, and an infinitesimal increase in 
temperature will result in a jump in reaction temperature from 
T-j. to Tg. At point I, which is defined as the "critical igni- 
tion condition" , the reaction moves from the kinetic regime 
to the diffusion regime of combustion. Temperature T^. is 
referred to as the ignition temperature. Although the des- 
cription of the process just presented is for a surface in 
which only heat transfer by convection occurs, the underlying 
ideas carry over for the total heat transfer at a point in the 
porous medium. The heat transfer can include contributions 
from conduction and radiation, as well as from convection. 

In his original work, Semenov [14] developed his theory 
to explain reaction behavior within homoegenous gas mixtures. 
Frank-Kamenetskii [12] later adapted this theory to describe 
the behavior associated with heterogeneous mixtures (i.e., a 
solid surface and air) . In their study of coal combustion. 
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Thomas, Stevenson, and Evans [15] state that a "temperature 
jump", as described by Frank-Kamenetskii , may or may not 
occur depending on the partial pressure of oxygen. This pres- 
ent investigation does not consider this aspect of the com- 
bustion theory. However, analyses were made using the 
present model to determine whether the reaction takes place 
in the kinetic regime or in the diffusion regime. Moreover, 
the Semenov model will be used to explain these results in 
Chapter V. 



D. ARRHENIUS LAW OF REACTION RATE 

For carbon reacting in air, Parker and Hottel [16] pro 

posed the following Arrhenius expression for the reaction 

2 

rate in units of kg-carbon/m -s. 



R 



= 9.55x10 



6 P °: 






(III .15) 



where P_ is the partial pressure of oxygen in atmosphere 
°2 

units, and is the universal gas constant (1.9 86 cal/gmole- 
K) . Equation III. 15 assumes a simple first order reaction 
for the combustion of carbon in air. Frank-Kamenetskii [12] 
has shown that the experimental data of Parker and Hottel is 
better correlated by a fractional order reaction, where the 
reaction order, n, is between 1/3 and 2/3. Replotting Parker's 
and Hottel' s data for n = 1/2, yields the reaction rate ex- 
pression used in the present analysis. 
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(III. 16) 



R = 2.065 xlQ 6 (R (j>) 1/2 exp( ~- 7 ?- ) 

c °2 R T 

u c 



2 

In expression III. 16 , R is in units of lbm-carbon/f t - 

c 

hr, R n is the gas constant for oxygen (48.29 f t-lbf/lbm-R) , 
u 2 

is in lbm/ft , R^ is 1.986 Btu/lbmole-R, and T c is in Rankines . 

The ideal gas law was used to replace the partial pressure, 

P n . Figure III. 3 and Table III.l show the results of plotting 
°2 

Parker's and Hottel ' s data for values of n equal to 1, 2/3, 

1/2, and 1/3. The values presented in Table III.l are in S.I. 
units . 

In order to determine the rate of heat generation and the 
rate of oxygen consumption, the chemical reaction for the 
combustion process must be considered. For nonvolatile com- 
bustion of carbon and oxygen, two reactions that describe the 
process are. 



c + | ° 2 -*■ CO (III. 17) 

C + 0 2 C0 2 (III. 18) 



The ratio of the mass rates of carbon monoxide to carbon 
dioxide produced increases with increasing temperature. 
Arthur [17] presents an expression for the rate ratio as a 
function of temperature (in Kelvins) . 



co/co 2 



2500 



,-6240 n 
exp (— ) 



(III. 19) 
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FIGURE III. 3 Results of plotting reaction rate 

data of Parker and Hottel [16] for 
several values of n. 
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R 



c 




gr-c 

cnrvs 



n 


a r i 


far- cal'] 


Icm-a+m-sed 


_ mole j 


i 

3 


1.19 x 10 3 


27500 


1 

2 


1.29 x 10 4 


31800 


2 

3 


2.1 x 10 5 


37100 


1 


9.5 5 x 10 6 


44000 



TABLE III.l Results of plotting reaction rate 

data of Parker and Hottel [16] for 
several values of n. 
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As a result of this dependence on temperature, the stoichio- 
metric ratio and the heat of reaction will also be functions 
of temperature. Defining the fraction of carbon monoxide 
being produced by, 



( CO/CO 2 ) 

F C0 = 1 + (co/co 2 ) 

and the fraction of carbon dioxide as, 

F = 1 
r C0 2 1 + (C0/C0 2 ) 

the heat of combustion is then expressed as. 



(III. 20) 



(III. 21) 



4H r - F co iH co + F C0 2 4 H co 2 (III - 22) 

Values for the heats of combustion, AH CQ and AH CQ , as func- 
tions of temperature were obtained from the JANAF (Joint Army, 
Navy, and Air Force) Tables [18] . The stoichiometric ratio 
(fuel to oxygen) of the overall reaction is. 



■R 



= f 



CO CO 



/(f 



CO. 



CO 



+ f 



CO CO. 



(Ill . 23) 



where f CQ is the stoichiometric ratio for reaction III. 17 and 
f^Q is the stoichiometric ratio for reaction III. 18. The 
rate of heat generation, R , and the rate of oxygen consumption 
can now be expressed by. 



34 



R 




(III. 24) 



g 



R 





(III. 25) 



It should be pointed out that Parker's and Hottel ' s work [16] 
and Arthur's work [17] were accomplished using specific types 
of carbon. Smoot and Pratt [19] and Frank-Kamenetskii [12] 
list tables and references for many other rate expressions, 
depending on the particular type of carbon used. The present 
model was developed for any rate expression of Arrhenius type 
for carbon consumption. For this investigation, we will re- 
main with Parker's and Hottel ' s expression III. 16 as modified 
by Frank-Kamenetskii. 

As previously stated, the particle diameter decreases as 
combustion progresses. The rate of decrease depends on the 
amount of carbon consumed at a point over time. Past observa- 
tions have shown that the effect is significant when the reac- 
tion is concentrated in a small region of the porous medium. 

To take this into account, expressions for the time rate of 
change of the diameter as a function of reaction rate were 
derived. The equation for cylindrical fibers is. 



d 



-2 R z D 2 /(tt p d) 
c c 



(III. 26) 



and for spherical particles 



d 



-2 R c z D 3 / (it p c d 2 ) 



(III. 27) 
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where p c is the bulk mass density of carbon. The diameter 
and the reaction rate are functions of both time and position. 

E. HEAT TRANSFER EQUATIONS FOR POROUS MEDIA 

In a previous investigation of porous media. Green and 
Perry [20] developed heat transfer equations for the solid 
and fluid phases. The model included the three basic mechan- 
isms, (1) physical movement of the fluid with its own heat 
content, (2) conduction of heat through the solid and fluid 
phases, and (3) convective heat transfer between the solid 
and the fluid phases. Radiation heat transfer was assumed 
to be negligible, and combustion was not a consideration. 

The differential equations for the solid and the fluid phases, 
obtained from energy considerations, 

3T 3 2 T 

- k s (1 -e>— r + hz(T w - T s» (III. 28) 

3 x 

3T 3T 3 2 T 

p c p-r— = -pp c U— r— — + k p ^ - hz(T - T ) (III. 29) 

w w^ 3t r w w 3x w^ £ X 2 w s 

and were solved numerically by Green and Perry using the 
finite difference method. In equations III. 28 and III. 29, 
c is specific heat, h is the internal convection heat trans- 
fer coefficient, k is thermal conductivity, and z is the wetted 
surface area per unit volume. The s and w subscripts refer 
to solid and fluid properties, respectively. Their model does 
not account for the temperature dependency of the properties. 
Other investigators such as Sundaresan [21] , Riaz [22] , 
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Mendelsohn [23], and Schneider [24] have presented a single 
equation approach to describe heat transfer in porous media, 




^-5 (III. 30) 
3x 



Equation III. 30 can be derived from equations III. 28 and 
III. 29 by assuming the solid and fluid phase temperatures are 
equal. Green and Perry proposed that the equal temperature 
assumption was valid for systems satisfying the following 
condition. 



where a is the thermal diffusivity of the fluid, 
w 

In the present investigation, the heat transfer equations 
are generalized to also include (1) radiation, (2) internal 
combustion, (3) temperature dependency of properties, and 
(4) compressibility effects of the air. 

An energy balance on the carbon gives the heat transfer 
equation. 



where k r is a pseudo thermal conductivity to account for 
radiation between particles. The derivation of equation 
III. 32 is presented in Appendix A. The effective conductivity, 
k g , of the porous solid was proposed by Russel [25] , 



X 




(III. 31) 




( 1 -P )p c c cTT (III * 32) 
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(III. 33) 



k 



e 



k 




where k and k are the bulk thermal conductivities of carbon 
c a 

and air, respectively, and p' = (1-p) . Russel's expression, 
based on an electrical resistance analogy, can be used for 
the full range of porosity, from 0 to 1. 

The pseudo thermal conductivity due to radiation, k^, is 
obtained as follows. Treating the idealized geometry of the 
porous medium as a series of closely spaced walls, the net 
radiation heat flux between two of the walls is. 



where e is the emissivity, and a is the Stef an-Boltzman con- 



stant. By expanding in a Taylor series about an 

analogy to Fourier's law of heat transfer by conduction yields 
the effective radiation conductivity as. 



Expression III. 35 is then used for k^ in equation III. 32. The 
details of the above derivation are presented in Appendix B. 
The application of expression III. 35 to porous media is also 
proposed by Rohsenow and Hartnett [26] . In a more rigorous 
development, Whitaker [27] presents an alternative approach 
for treating radiation heat transfer in porous media. 



oe(T* - T^) / (2-e) 



(III. 34) 
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4 a e 6 T / (2-e) 
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Following the experimental work of Yoshida, Ramaswami, and 
Hougen [28] , the internal convection heat transfer coefficient 
is given by the empirical correlation, 

h = 0.91 Re ,-0 * 51 U> c G (c y/k ) ~ 2/3 ] (III. 36) 

a a a fm 

where ip is equal to .91 for cylindrical fibers, and equal to 
1 for spherical particles, G is a pseudo mass velocity given 
by pp u, and c is the specific heat of air at constant pres- 
sure. The fm subscript refers to the properties evaluated 
at film temperature. Re' is a pseudo Reynolds number defined 
by. 



Re ' = G/(z v \p) 



(III. 37) 



The air properties vary with temperature, and h varies with 
temperature. Finally, the reaction rate for the heat genera- 
tion term in equation III. 32 is given by expression III. 24. 

An energy balance on the air within the porous medium 
provides the second heat transfer equation as, 



9T a 9T 
^ (pk alF> ‘ Pf:> a C a u 3x 



a + hz (T - T ) + U 3 <P P) . 
c a 3x 



3T 

P p a C aTF 



(III. 38) 



The density of air, p , is approximated in the model by the 

3 

ideal gas law. The term, 3(pP)/3x, is due to the compressi- 
bility of the air. The derivation of equation III. 38 is 
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presented in Appendix A. All properties in equation III. 38 
vary with temperature . The properties of standard air were 
used in the analysis. During combustion, oxygen in the air 
is replaced by carbon monoxide and carbon dioxide. Calcula- 
tions showed that the gas mixture would have thermophysical 
properties slightly different from those of standard air. 

In a "worst case" situation, the average difference between 
the thermophysical properties of air and a mixture of .79 ^ 
and .21 CC >2 is 7 percent. This mixture assumes that the oxy- 
gen is totally consumed. Viscosity had the largest difference 
of 13 percent, and specific heat the smallest at 1.5 percent. 

The properties of a .79 ^ and a .21 CO mixture are approxi- 
mately those of air (different by less than 2 percent) . At 
typically observed temperatures (greater than 1200 deg-F) , the 
combustion gas will be mostly carbon monoxide, and the error 
introduced by assuming the properties of standard air is mini- 
mum. A small difference in the properties is acceptable since 
an additional molecule mass transfer equation for either carbon 
monoxide or carbon dioxide would be necessary to calculate 
the effective properties of the gas mixture. The methods used 
to obtain the thermophysical properties of the gas mixtures 
above were those of Mason and Saxena [29] , and Reynolds and 
Perkins [30] . The polynomial expressions used to calculate 
the thermophysical properties of air are derived in Appendix B. 

F. OXYGEN DIFFUSION EQUATION FOR POROUS MEDIA 

As a result of combustion, the heat transfer equations 
include four response variables, the carbon and air temperatures. 
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T and T , the concentration of oxygen, <f>, and the total 

C cl 

internal pressure, P. The fourth field equation necessary 
to complete the system is obtained from species conservation 
considerations. The oxygen molecule transport mechanisms 
included in the model are (1) molecular diffusion resulting 
from concentration gradients (Fick's Law), (2) convective 
mass flow, and (3) oxygen consumption due to combustion. 
Diffusion resulting from pressure and temperature gradients 
was considered negligible. An example cited by Bird, Stewart, 
and Lightfoot [21] where pressure diffusion is important is 
the centrifuge separator; and for thermal diffusion, the 
Clusius-Dickel column, where very steep temperature gradients 
are used to separate complex mixtures of organic molecules. 
Thus, the oxygen molecule transfer equation for a porous 
medium is. 



3x (pI? e 3x^ “ 3x (pU< ^ ) " R ° 2 2 P 3t (III. 39) 

The derivation of equation III. 39 is presented in Appendix A. 
The second order term in equation III. 39 results from Fick's 
law of diffusion. The effective diffusivity of the porous 
medium, V , is a function of pressure, temperature, and pore 
geometry of the medium. Scheidegger [9] presents several 
models for V which have been proposed by a number of inves- 
tigators. Briefly, there are two limiting cases of diffusion, 
(1) diffusion associated with the collision between gas mole- 
cules, and (2) diffusion associated with the collision of 
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gas molecules with the pore walls of the medium. This second 
phenomenon is discussed by Bennett and Myers [32] , and is 
referred to as Knudsen diffusion. The former case occurs 
when the mean free path of the molecules is smaller than the 
pore diameter. The expression for mean free path, gj, given 
by Treybal [33] is, 

oj = ^^-(R u T/2 u g c M) 1/2 (III. 40) 

where g^ is the gravitational constant and M is the molecular 
weight. From equation III. 40, a value of co = 6/100 is obtained 
for oxygen for the geometry of typical porous media encoun- 
tered in the analysis. Thus, collision between molecules is 
the governing diffusion mechanism. A semi-empirical expres- 
sion, proposed by Gilliland [34], is used to obtain the diffu- 
sion coefficient of oxygen into air. Gilliland's expression 
is given by, 

V = 435.7 T 2 ^ 3 (M _1 + m" 1 )/ [P (V 1 / 3 + vi /3 ) ] (III. 41) 

cl Si u 2 3 . u 2 

2 

where V is in units of cm /sec, P is the total pressure in 

Pa, V and V_ are the molecular volumes of air and oxygen, 
a o 2 

respectively. M and M n are the molecular weights of air and 

a ° 2 

oxygen. The values of V and V were obtained from Holman 

a C> 2 

[35] as 29.9 and 7.4, respectively. The effective diffusivity 
proposed by Denbigh and Turner [36] for a porous medium is. 
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V = V/t (III. 42) 

Expression III. 42 accounts for the tortuous path the oxygen 
molecules follow through the porous medium. Lastly, the oxy- 
gen consumption term in equation III. 39 is given by expression 
III. 25. 

G. THE SURFACE RECESSION PROBLEM 

The previous system of equations describes the combustion 
process occurring within the porous medium. As combustion 
progresses, the oxygen concentration goes to zero, and as a 
result, the reaction moves to the air inlet surface of the 
porous medium. Thereafter, the carbon consumption takes place 
at the surface. During the surface recession phase, the thick 
ness of the porous medium decreases with time, i.e., surface 
recession. This suggests partitioning the problem into two 
parts, (1) initial combustion within the porous medium, and 
(2) combustion at the air inlet surface. The second phase 
of the problem is formulated as a moving boundary problem. 

The governing equations for the moving boundary problem are 
those developed previously with the exception of the oxygen 
molecule mass transfer equation. With combustion occurring 
at the surface, oxygen is totally consumed in a shallow region 
near the surface. Experimentation by Koizumi [37] has shown 
this penetration of the oxygen to be 1 to 2 particle diameters 
As a result, the oxygen molecule mass transfer equation is 
not needed for this phase of the problem, and the heat genera- 
tion can be treated as a planar source. 
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The carbon heat transfer equation with the modified heat 



generation term is, 



l- 2 §-[ u-p) < k e+ k r ) 



3T 



] -hz (T -T ) +R' z5 A (ti = 0) 
c a g 



c 



r 3n 




(III. 43) 



where n is equal to x/L. The first term on the right side of 
equation III. 43 arises from the x coordinate becoming a func- 
tion of time during surface recession. A similar term appears 
in the air temperature and combined Darcy's law-continuity 
equations as well. Transformation of the field equations from 
a fixed coordinate to a moving coordinate system is presented 
in Appendix A. As was observed (see Figure V.2) , the tempera- 
ture gradient through the porous medium during surface reces- 
sion was small. As a result, the additional term may be 
omitted without affecting the results. The moving boundary 
problem formulation is based on Crank’s [38] extension of a 
method proposed by Murray and Landis [39] . The thickness, L, 
is updated continuously by evaluating its time rate of change 
given by. 



Expression III. 44 is obtained from a mass balance on the 
carbon (i.e., the time rate of change in carbon is equal to 



L 



-(puf ) <p /[(1-p) 

x=0 x=0 




(III. 44) 
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the carbon consumption) . Noting that all the oxygen entering 

the plate is consumed in a small region near x = 0, the carbon 

consumption can be represented in terms of the oxygen flow 

and the stoichiometric ratio (i.e., R' = -puLif A where A 

c r R 00 

• 

is area) . The time rate of change for carbon is (l-p)p c AL. 

The planar heat source, R^ , in equation III. 43, is determined 
by the amount of oxygen entering the porous medium at x = 0. 
This is converted to heat generation by using the stoichio- 
metric ratio, f.,, and the heat of combustion, AH . Thus in 
equation III. 42, the planar heat source becomes puf^H^^. 
Ozisik [40] presents a general discussion and references for 
the moving heat source problem. 

Transition from the internal combustion problem to the 
surface recession problem occurs when the porosity at x = 0 
nears a value of 1. The computer program, implementing the 
analysis, was written so that the change from an internal 
combustion problem to a surface recession problem occurs 
automatically . 

H. BOUNDARY CONDITIONS 

Three sets of boundary conditions may be used for equations 
III. 32, III. 38, and 11.39. Each set of boundary conditions 
approximates a physical situation (e.g., thermal flow reactor 
or Fontenot's experiments [1]). The boundary conditions 
were the nearest approximations that could be made and remain 
within the limitations of a one dimensional model. The first 
and second set of boundary conditions are typical of thermal 
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flow reactors with and without radiation from the boundary 



surfaces. These are 
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and with radiation. 
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x = 0 (III. 45) 
x = L (III. 46) 
x = 0 (III. 47) 
x = L (III. 48) 
x = 0 (III. 49) 
x = L (III. 50) 

x = 0 (III. 51) 
x = L (III. 52) 
x = 0 (III. 53) 
x = L (III. 54) 
x = 0 (III. 55) 
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0 



x = L 



(III. 56) 



34 > 

3x 

Expressions III. 45 and III. 46 represent insulated boundaries 
on the porous solid (i.e., no heat loss from the porous solid 
to the environment). Expressions III. 51 and III. 52 provide 
for radiation heat transfer from the porous solid boundaries 
to the environment. The boundary conditions for the air 
temperature and oxygen concentration, expressions III. 47 
through III. 49 and III. 52 through III. 55, are the Danckwerts ' 
boundary conditions for flow reactors. A convincing discussion 
of the Danckwerts' boundary conditions applied to mass diffu- 
sion equations for porous media is given by Bischoff [41]. 

A brief summary of the discussion was presented by Vatikiotis 
[2]. An analogy to Bischoff's approach for the air heat trans- 
fer equation is presented in Appendix B. 

The third set of boundary conditions approximates the situa- 
tion for Fontenot's experiments [1]. In the experiments, plate 
specimens were mounted into the wall of a wind tunnel. The 
interior surface of the plate (i.e., at x = L) was exposed to 
the wind tunnel air flow, and the external plate surface (i.e., 
at x = 0) was exposed to the outside environment. The expres- 
sion for calculating the pressure differential across the 
porous medium as a result of the external flow over the surface 
at x = L is formulated in Appendix B. It must be emphasized 
that the purpose for modelling Fontenot's experiments was to 
analyze the combustion behavior within the plate specimens 
using the one dimensional model. However, the external flow 
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and plate surfaces involved convection heat transfer requiring 
two dimensional analysis. Noting this, the following approach 
was adopted to provide for convection heat transfer at the 
boundaries of the one dimensional model in a qualitative 
sense. The boundary conditions are. 



(k e +k r ) TF = h 1 (T c -T oo )+a £ (T c -T o3 ) 
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x = 0 (III. 57) 

x = L (III. 58) 
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x = 0 (III. 59) 




0 



x = L (III. 60) 



(p = . (p m x = 0 (III. 61) 

|^=0 x = L (III. 62) 

3x 



As stated above, the values of h^ and h ^ in the convection 
terms of expressions III. 57 and III. 58 are difficult to inter- 
pret for a one dimensional model. The difficulty arises from 
considering the boundary layers on the external surfaces of 
the porous medium (e.g., considering the porous medium as a 
flat porous plate, the boundary layer resulting from natural 
convection on the x = 0 surface increases in thickness in the 
vertical direction, and for forced convection at the x = L 
surface, the boundary layer thickness increases in the flow 
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direction) . For this investigation, was approximated 
by. 



k up 

h i - < \ a > 0 

x=0 


(III. 63) 


This expression is based on the analytical results 


presented 



by Merkin [42]. Expression III. 63 is obtained by considering 
natural convection heat transfer from a vertical flat plate 
with "suction". From Merkin ' s results, the boundary layer 
thickness becomes constant at some distance in the vertical 
direction, and thus, the convection heat transfer may also be 
considered constant. The magnitude of the convection heat 
transfer coefficient, h 2 , varies in a direction parallel to 
the external flow, U at the x = L surface. In addition, h„ 

oo z 

depends on the efflux of the gas at the surface. To simplify 



the analysis, h 2 was approximated by the relations 
flat plate given by Holman [35] as. 


for a smooth 


h 2 = .664 -p Pr 1/3 Re 1/2 


(III .64) 


for laminar flow, where L is the distance from the 
layer initiation, and. 


boundary 


h 2 = Pr 1/3 ( .037 Re* 8 - 850) 


(III. 65) 


for turbulent flow, where Pr is the Prandtl number. 


. The 


Reynolds number here is based on the external flow 


velocity 
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over the surface at x = L. Kays [43] provides an alternate 
scheme for treating boundary layers with "blowing and suction". 
The air temperature and oxygen concentration boundary condi- 
tions, III. 59 through III. 62 were also selected because of 
the difficulties in treating external surface boundary layers 
in a one dimensional model. Alternatives to the essential 
boundary conditions of the air and oxygen are the Danckwerts ' 
conditions. However, using the Danckwerts' conditions with 
convection boundary conditions poses the problem of approxi- 
mating a reasonable value for u at x = 0 and x = L (i.e. , 
boundary layers affect the component of velocity, u, normal 
to an external surface). To overcome these difficulties, it 
seemed prudent not to have pore velocity and the convection 
heat transfer coefficients, h^ and , appear in the same set 
of boundary conditions. However, the difficulty still re- 
mains of specifying a distance from the origin of the boundary 
layer for expressions III. 64 and III. 65. In the initial work 
by Vatikiotis [2] , an arbitrary distance of 1 foot was used 
to compare results. 

I. INITIAL CONDITIONS 

To initiate combustion, the porous medium must be brought 
to a temperature at which the heat generation rate is suffi- 
ciently high. Experimentally, there are many techniques for 
doing this. However, it is difficult to measure the actual 
temperature and oxygen concentration during the experiments. 
This difficulty is carried over to identifying a reasonable 
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set of initial conditions. Noting this, the model was de 
veloped so a problem can be started at ambient conditions 



with a heat flux applied to the carbon at a boundary 
heat flux is mathematically treated as, 

3T 

(1 -p> <k e + V = -%t 

This approach of starting the problem simplifies the 
determining a reasonable set of initial conditions, 
tion, any arbitrary set of initial conditions can be 



The 

(III. 66) 

task of 
In addi- 
specif ied. 
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IV. IMPLEMENTATION OF NUMERICAL METHODS 



A. SOLUTION ALGORITHM 

The field equations and auxiliary equations which define 
the problem were presented in the previous sections. The 
following is the scheme used to solve for the dependent 
variables T , T , 4>, u, P, L, and p : 

C a 3. 

1. Starting with the Dupuit-Forcheimer expression III. 10, 
the filter velocity, Q, is represented as a function 
of the pore velocity, u. 

2. Substitute for Q from the previous step into Darcy's 
law, equation A.l. 

3. The expression formed in step 2 is solved for the pore 
velocity, u, as a function of dP/dx. This procedure 
yields equation A. 2. 

4. Equation A. 2 is substituted into the continuity equation 
A. 3 resulting in the 2nd order partial differential 
equation A. 4 with the pressure, P, and the air density, 

p , as the dependent variables. 

a 

5. Expanding equation A. 4 yields equation A. 5. Equation 
A. 5 is one of 4 equations which were integrated to 
produce the problem solution. By taking the air den- 
sity, p , as a constant over a time step of integration, 
equation A. 5 becomes an ordinary differential equation 
with respect to x. For equation A. 5, the shooting method 
(with Euler's formula for the integration algorithm) 
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was used to solve for the pressure, P(x) , and the 
pressure gradient, dP(x)/dx at each time step. The 
air density, p , was updated at each time step using 
the ideal gas law, equation B.41. The rate of change 
of the air density was neglected during the first two 
time steps of integration. Reasons for neglecting the 
time rate of change of the air density are given in 
Appendix C.2. 

6. Having solved equation A. 5 for the pressure gradient, 
dP/dx, the pore velocity can be obtained from equation 
A. 2. 

7. In solving the air heat transfer equation III. 38 and 
the oxygen molecule transfer equation III. 39, the pore 
velocity, u, and the air density, p , are obtained from 
equations A. 2 and B.41, respectively. Equation III. 3 8 
and III. 39, along with the carbon particle equation 
III. 32, were transformed to ordinary differential equa- 
tions in time by a Galerkin FEM formulation. 

8. The ordinary differential equations from step 7 are 
given by equations C.17 through C.19. The time inte- 
gration was performed by a version of Gear's method 
developed by Franke [44] . 

9. All thermophysical properties (i.e., k e , k^, k , c a , 

c , y) being functions of temperature are continuously 
updated during the transient analysis. 

10. As carbon is consumed during combustion, the porous 

medium properties (i.e., p, m, 5, d) change with time 
and are updated continuously as well. 
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12. No attempt was made to match the orders of convergence 
of the integration algorithms. 

Implementing the above solution algorithm resulted in a 
computer program for the solution of the problem. Several 
options were incorporated in the computer program allowing 
the user flexibility in the analysis. Boundary conditions 
simulating fixed air inlet temperature, and convection and 
radiation heat transfer from the exterior surfaces can be se- 
lected. The problem can be initiated in two ways; an arbitrary 
set of initial conditions can be specified, or a heat flux, 
shown by expression III. 66, can be specified with the porous 
medium at ambient conditions. In addition, the porosity at 
x = 0 causing transition to the surface recession phase can 
be specified. Other options permit investigating materials 
other than carbon (e.g., boron), studying the effects of non- 
constant distributions of porous media properties. 

The computer program includes a main program and 24 sub- 
routines (7 of the subroutines comprise Franke's [44] inte- 
gration routine) . The function and description of each of 
the subroutines are found in the computer program listing. 
Figure IV. 1 shows the flow of the program, and briefly des- 
cribes the process at each step. Each block represents a 
subroutine. The configuration of the blocks (i.e., blocks 
within blocks) in Figure IV. 1 show the six levels of the 
computer program. In addition, other processes associated 
with the integration routine occur inside the block shown by 
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Main program. 



FALSE- 



Read in parameters, initialize -flags and 
counters; Calculate initial properties 
and pore velocity. 

f 



Construct FEM nodal poi nt/el ement 
correspondence array. 



Print current response variables 
and properties. 



Check -for transition 


TRUE 


Reset -flags, counters delete 


to surface recession. 




<t> (called once). 


T 

FALSE 








t 






L_ 



Call integration routine -for T , T 3 , (p , ( (p deleted 

during surface recession). 


^Construct O.D.E.’s. 
1 


1 

1 


Calculate properties. 


1 

1 




Test -for oxygen instabiltiy. 


1 

| 




I-f surface recession, update thickness, L. 


1 

1 

1 




Calculate pore velocity and pressure gradient. 


1 

1 

1 




I-f surface recession, trans-form reaction rate 
expression to surface -flux. 


1 

1 




Set -flag -for sur-face recession i -f p(n=l) = pm ax. 


1 






1 

1 

1 


Zero-out system matrices. 


1 

1 

1 


Generate FEM operators. 


1 

1 

| 


Form mass and system matrices. 


1 

1 

1 


Adjust matrices -for boundary conditions. 


1 

1 




Adjust matrices -for oxygen instability. 


|Per-form integration. 



-Test -flags -for program termination. 



1 

TRUE 

!_ 



Print last values o-f response variables 
and properties. 



T 

STOP 



FIGURE IY.l - Flow 
of computer program. 
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the dashed line. These processes are discussed in detail in 
Franke's report. 

B. STIFFNESS CONSIDERATIONS 

Initially, the ordinary differential equations resulting 
from the finite element formulation were integrated explicitly 

using a sixth order Runge-Kutta method (IMSL subroutine DVERK) . 

-4 

The largest obtainable time step was approximately 10 
seconds. A larger time step caused the integration to become 
unstable. Gear [45] described this behavior as typical of 
"stiff" differential equations, and suggested that implicit 
integration methods be used to improve the efficiency of the 
integration. Gear's method, presented by Brown and Gear [46] 
and modified by Franke [44], was adopted as the integration 
method. As a result, time steps of several seconds were ob- 
tained for some problems (i.e., for surface recession and when 
extinction occurred) . In addition, moving the exponential 
reaction rate from the excitation vector to the stiffness 
matrix (discussed in Appendix C.l) improved the integration 
efficiency. The improved efficiency was attributed to a better 
approximation of the Jacobian matrix in Gear's method. 

Bui [47] and Burka [48] discuss the difficulties of solving 
differential equations exhibiting "stiffness" and propose 
alternatives to Gear's method. For a system exhibiting "stiff- 
ness", small time steps are needed for the rapidly responding 
terms while the integration must be continued for a long period 
to account for the slowly responding terms . A measure of 
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"stiffness" for a system of differential equations is the 
ratio of the real part of the maximum eigenvalue of the Jacob- 
ian matrix to the real part of the minimum eigenvalue. This 
is called the "local stiffness ratio", and is represented 
by. 

Max | A . | 

n __ 1-1 , ... , n 

Min 

i—l , ... , n 

where n is the number of equations. A system of equations can 
be considered "stiff" if S has a value greater than 10, and 
the real parts of the eigenvalues are negative. 

Lastly, most numerical integration methods are bound by 
a maximum time step for which the solution remains stable. 

This is defined by. 



x i 






| At • A ± | < K 



(IV. 2) 



where At is the time step, A^ is the real part of the maximum 
eigenvalue of the Jacogian matrix, and K is a constant depend- 
ing on the integration method (usually 1 < K < 10) . For 
example, Bui [47] states that Euler's method requires 
j A t A ^ | < 1. Therefore, for a maximum Jacobian matrix eigen- 
value (representing fastest response) of 50 (seconds) the 
largest timestep for the integration is .02 seconds. Exceed- 
ing this value results in the integration becoming unstable. 
This explains the maximum timestep of 10~^ seconds for the 
sixth order Runge-Kutta method as previously discussed. 
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C. TREATMENT OF NUMERICAL DIFFICULTIES ASSOCIATED WITH 

OXYGEN CONCENTRATION 

In problems of combustion or kinetics, where dependent 
variables go to zero, numerical difficulties arise. In his 
discussion, Frank-Kamenetskii [12] points out that for frac- 
tional order reactions (expression III. 16), both the oxygen 
concentration and concentration gradient can become equal to 
zero within the porous medium. This requirement is difficult 
to satisfy with approximate solution methods. In the initial 
effort (Vatikiotis [2]), rim time errors resulted when the 
oxygen concentration became slightly negative. This was 
caused by attempting to take the logarithm of a negative con- 
centration in the fractional order term in expression III. 16. 

One approach for correcting this problem was to use the abso- 
lute value of the oxygen concentration in the reaction rate 
term. This resulted in the concentration becoming increasingly 
negative. Another approach was to check the values of oxygen 
concentration at each integration interval and setting the 
concentration to zero for the negative values. At the next 
integration, the reaction rate was zero at the locations of 
zero concentration. Without consumption the oxygen concentra- 
tion at the nodal point previously set to zero becomes positive 
and the oscillation repeats itself. It was observed that de- 
creasing the integration time step and the length of the ele- 
ments (i.e., increasing the number of nodal points) minimized 
the oscillations. Since CPU time increased for a given problem, 
the approach was not totally satisfactory. However, this 
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method was adopted for obtaining the results presented by 
Vatikiotis [2] . 

An attempt to resolve the above numerical difficulty for 
the present model was by implementing the Moving Finite Element 
method of Gelinas, Doss, and Miller [49]. In the method, 
the nodal points defining the elements converge to regions 
where they are needed to minimize the error in the numerical 
solution. As a result, the size of the elements in a region 
of large reaction rate become small, and move with this region 
in time. The motivation for implementing the Moving Finite 
Element method was Frank-Kamenetskii ' s [12] discussion of 
fractional order reactions as stated above. Also, it was ob- 
served that the oxygen concentration approached zero in a 
small interval within the porous medium. This behavior resulted 
in steep concentration gradients. It was thought that a num- 
ber of small elements in this region would improve the accuracy 
of the solution. The results obtained from the Moving Finite 
Elements showed the method was not effective in resolving the 
numerical difficulties discussed above. The method was subse- 
quently abandoned because of the additional expense associated 
with solving the differential equation governing the nodal 
point locations. 

The method adopted in the present model for minimizing 
the difficulties with the oxygen concentration is as follows. 

The time derivative of the oxygen concentration at a nodal 
point is set to zero when the concentration reaches a speci- 
fied positive minimum value (less than 10 ^ lbm/ft^) . At 
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the next integration interval, the time derivative remains 

zero if the concentration at the adjacent nodal point in the 

negative direction (i.e., upstream) has decreased. If the 

concentration at the adjacent nodal point has increased, then 

the time derivative is released. This method, in conjunction 

with Franke's [44] integration routine, appears to act in 

an iterative manner continually improving the solution. The 

values of concentration in the regions where the oxygen has 

been totally consumed are typically on the order of 10 ^ 

3 

lbm/ft . Although they do not give the details, Scheisser 
and Stein [50] have used the above method in their work of 
simulating coal conversion. 

D. AN ALTERNATE SOLUTION STRATEGY 

The experience and difficulties of implementing the numeri- 
cal methods, as discussed above, have provided insight for an 
alternate strategy for solving the combustion and heat trans- 
fer problem. The alternate strategy consists of formulating 
the problem with two moving boundaries. One moving boundary, 
as in the present formulation, would account for surface 
recession (i.e., the porous medium thickness changing with 
time) . The second moving boundary would account for the 
oxygen concentration going to zero within the porous medium. 
This second boundary would be located where the concentration 
becomes zero, between x = 0 and the first moving boundary. 

The advantages of the alternate solution strategy are (1) 
eliminating the numerical difficulties associated with the 
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oxygen concentration (discussed in the previous section) , 

(2) relaxing the planar heat source assumption for the surface 
recession phase (discussed in Section III.G), and (3) possibly 
improving the numerical "stiffness" characteristics of the 
integration (discussed in Section IV. B) . 
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V . RESULTS AND OBSERVATIONS 



A. BASIS FOR THE RESULTS 

The results presented in the following sections are based 
on the idealized porous medium described in Section III. A and 
Figure III.l. For the idealized porous medium, porosity and 
permeability were explicit functions of pore diameter and unit 
cell thickness. Other porous media will have different geo- 
metric relationships for describing the properties. In addi- 
tion, for naturally occurring porous media, experimental methods 
would be employed to determine the geometric properties. 

Noting this, the results of this investigation are intended 
to describe, in a qualitative sense, the behavior of a carbon 
porous medium during combustion. As will be shown, the tempera- 
tures observed during combustion are highly dependent on the 
geometric properties and the pore velocity. Although, the 
temperatures are typical of those reported for carbon combustion 
experiments, it would not be reasonable to perform a quanti- 
tative comparison between the results obtained by the model 
and those obtained experimentally unless the porous medium 
properties were the same for both. This point will be demon- 
strated in the next section. 

B. COMPARISON OF COMBUSTION MODEL RESULTS TO EXPERIMENTAL 

RESULTS 

The behavior of the temperature and oxygen concentration 
of the combustion and heat transfer model is in fair agreement 
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with the behavior reported in the experimental investigations 
of Vulis [13] , Koizumi [37] , Fontenot [51] , Spalding [52] , and 
Kolodstev [53] . A quantitative comparison of the results can 
not be made for two reasons. First, there have been few 
experimental investigations which present data for transient 
behavior of combustion in porous media. Most of the analyses 
present quasi-steady data (quasi in the sense that for some 
period of time, there is little change in the response varia- 
bles) . When transient behavior is reported, it is usually 
presented in a narrative form. Secondly, the reported des- 
criptions of the porous media used in the experiments are not 
sufficient. The properties of the carbon used (e.g., thermal 
conductivity, density, Arrhenius rate law, etc.) and the 
porosity and permeability relationships are needed to accurately 
simulate the experiment. The profiles of the quasi-steady 
temperatures obtained by the combustion model and those ob- 
tained by Kolodstev 1 s experiments will be shown together. The 
purpose of this is to demonstrate that although Kolodstev 's 
porous medium could only be roughly approximated, the tempera- 
ture profile obtained by the combustion model was similar in 
shape and magnitude to that obtained experimentally. 

In the experimental investigations of Vulis [13] , Fontenot 
[51], and Spalding [52] , it was observed that for a given 
temperature, sustained combustion occurred for a specific range 
of velocities. Velocities outside the range, either greater 
or lower, would cause the combustion to extinguish and the 
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porous medium to cool to ambient temperature. This phenomenon 
was also observed for the combustion model. The explanation 
of this behavior, discussed in Section V.G, is summarized as 
follows. As pore velocity increases, convection heat transfer 
and oxygen supply for combustion both increase. Similarly, 
a decrease in pore velocity will cause a decrease in convec- 
tion heat transfer and oxygen supply. For a given set of 
conditions, pore velocities above the range will result in 
the heat transfer or heat loss overcoming the heat generation. 
Also, for pore velocities below the range, the oxygen supply 
for combustion is decreased to the extent that the heat genera- 
tion is overcome by the heat loss. In other words, the higher 
pore velocities tend to "blowout" the reaction (e.g., blowing 
out a candle) , and the lower pore velocities tend to "starve" 
the reaction (e.g., placing a container over a candle) . 

Important at lower pore velocities are the effects of heat 
transfer from the external surface of the porous solid (i.e., 
at the boundaries) . The experiments performed by Vulis [13] , 
Fontenot [51] , and Spalding [52] allowed heat transfer from 
the surfaces by either convection or radiation. The effects 
of heat transfer at the boundaries on combustion are discussed 
in Section V.G. 

Kolodstev [53] investigated combustion gas dynamics in a 
porous medium comprised of spherically shaped carbon particles. 
This is the same type of carbon used by Parker and Hottel [16] . 
In his results, Kolodstev presents the quasi-steady air tempera- 
ture as a function of the depth of the porous medium (particle 
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temperature was not reported). Figure V.l shows Kolodstev's 
results plotted with the results obtained by the combustion 
model. Kolodstev's experiment could not be simulated accu- 
rately. As stated previously, the description of Kolodstev's 
porous medium was not sufficient. The description of the 
porous medium and the ambient conditions used to simulate the 
experiment is shown in Table V.l. Also shown are the known 
parameters given by Kolodstev. In Figure V.l, the temperature 
for the first inch of the 7.5 inches of the porous medium are 
shown (temperatures at greater depths were not reported) . As 
can be seen, there is fair agreement in the shapes of the 
profiles. In addition, the magnitude of the temperatures pro- 
duced by the combustion model are similar to those obtained 
in the experiment. 

The combustion model results showed that once the reaction 
moved to the air inlet surface of the porous medium, the oxy- 
gen concentration penetrated the porous medium by 1 or 2 
particle diameters. This observation was also reported in 
the experimental results of Koizumi [37] and Kolodstev [53] . 
The response of the oxygen concentration, in general, will be 
discussed in the next section. 

C. EXAMPLE RESULTS 

The intent here is simply to demonstrate the operation of 
the computer program. The behavior of a single porous medium 
subjected to several initial conditions is presented. The 
description of the porous medium and the ambient conditions 
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FIGURE V . 1 Air temperature from combustion model 

shown with that obtained from 
Kolodstev's [53] experiment. 



TABLE V. 1 Geometry o-f porous medium and ambient conditions 

■for simulating Kolodstev’s C533 experiment. 



• Particle Shape: 


spher i cal 


• Particle Diameter (d): 


.126 in 


Unit Cell Thickness ( D ) : 


.126 in 


• Spatial Thickness o-f Porous Medium ( L. ): 


7.5 in 


Porosi ty ( p ) 


. 476 


Per meab i 1 i t y ( m ) 


_c 2 

.103x10° -ft 


Bulk Thermal Conductivity o-f Carbon ( k c ) : 


86.0 Btu/-f t-hr-F 


Bulk Speci-fic Heat o-f Carbon ( C c ) : 


.231 Btu/lbm-F 


Bulk Density o-f Carbon (p ): 


70.3 lbrn/Tt 3 


Thermal Emissivity o-f Particles ( € ): 


. 9 


Ambient Temperature ( *J^) : 


80.0 deg-F 


• Ambient Pressure 


14.7 psi 


• Pressure Di -f -f erenti al (aP) : 


-.29 psi 


•Ambient Oxygen Concentration (C^) : 


.0172 1 bm/-f t 



• Known parameter -for Kolodstev's experiment. 
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used for the following problems are shown in Table V.2. The 
boundary conditions for the examples in this section were 
insulated boundaries on the porous solid temperature (ex- 
pressions III. 45 and III. 46), and the Dankwerts conditions on 
the air temperature and oxygen concentration (expressions 
III. 47 through III. 50). 

As the first example (referred to as example A) , a heat 
4 2 

flux of 3.0 x 10 Btu/ft -hr was applied at x = 0 for 15 seconds. 
The porous medium was initially at a constant temperature of 
80 degrees Fahrenheit and a constant oxygen concentration of 

3 

.0172 lbm/ft (ambient conditions) . The results of this prob- 
lem are shown in Figures V.2-V.4. Figure V.2 shows the tem- 
perature increasing after the heat flux was removed. At 25 
seconds, the porosity at x = 0 reached .95 and the problem 
transitioned to the surface recession phase. The values of 
thickness for the times shown in Figure V.2 are .25 inches to 
25 seconds, .226 inches at 72 seconds, .133 inches at 357 
seconds, and .0238 inches at 517 seconds. The computer pro- 
gram was written such that the problem terminates when the 
thickness becomes 10 percent or less of the initial value. 

The air temperature shown by the dashed line in Figure V.2 
follows the carbon temperature with the exception of the air 
inlet region at x = 0 . Figure V.3 shows the oxygen behavior 
for this problem. At 25 seconds, the oxygen concentration at 

-4 3 

x = 0 was 3.6 x 10 lbm/ft and the penetration of the oxygen 
was to x/L = .01. This supports the assumption of using a 
planar source for the heat generation and deleting the oxygen 
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TABLE V.2 Geometry of porous medium and ambient conditions 

■for the example problems in Section V.C. 



Particle Shape: 


spheri cal 


Particle Diameter (d): 


. 005 i n 


Unit Cell Thickness ( D>: 


. 005 i n 


Spatial Thickness of Porous Medium ( L>: 


.25 in 


Porosity ( p ) 


. 476 


Permeability ( rn ) 


-9 2 

.162x10^ ft 


Bulk Thermal Conductivity of Carbon ( k c ) : 


86.0 Btu/f t-hr-F 


Bulk Specific Heat of Carbon ( C c ) : 


.231 Btu/lbm-F 


Bulk Density of Carbon ( p c ) : 


70.3 lbm/ft^ 


Thermal Emissivity of Particles (€ >: 


. 9 


Ambient Temperature ( : 


80.0 deg-F 


Ambient Pressure (PL): 

CO 


14.7 psi 


Pressure Differential (Ap) : 

Ambient Oxygen Concentration < <i) : 


-.35 psi 
.0172 lbm/ft^ 
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FIGURE V. 3 Oxygen response during sustained combustion. 
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FIGURE V.4 Change in thickness of porous medium during 

surface recession phase. 



concentration equation during the surface recession phase. 

Figure V.4 shows the decrease in the thickness of the porous 
medium with time. The nondimensional parameter, X, given by 
expression III. 31, had an average value of .63 at the start of 
the problem. This value of X is above the minimum value of 
.342 proposed by Green and Perry [20] when assuming equal 
carbon and air temperatures. This suggests that a single heat 
transfer equation would have been sufficient to describe the 
temperature response for this problem. However, this is not 
valid for all cases. For example, when d, D, and L are 
doubled, that is, .01 inches, .01 inches, and .5 inches, 
respectively, X equals .22, and for d, D, and L equal to .025 
inches, .025 inches and 1.25 inches, respectively, X equals 
.056. Hence, for the general case, the carbon and the air 
temperatures should be treated as separate response variables. 

As a second example (referred to as example B) , a heat 
4 2 

flux of 3.0 x 10 Btu/ft -hr was applied at x = 0 for 14 seconds. 
The porous medium was initially at the same ambient tempera- 
ture and oxygen concentration as in the previous example. 

Figures V.5 and V.6 show the results of this problem. The 
porous medium cooled to ambient conditions in approximately 
41 seconds after the heat flux was removed. Applying the 
heat flux for 14 seconds was not sufficient to produce the 
conditions needed to sustain combustion. However, increasing 
the duration of the heat flux for an additional second results 
in sustained combustion. Figures V.5 and V.6 show the results 
starting at 12 seconds. The behavior of the system prior to 
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FIGURE V.5 Response of temperature when heat flux at x 

insufficient to produce sustained combustion 
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12 seconds is the same as in Figures V.2 and V.3 of the pre- 
vious problem. 

As the last example (referred to as example C) , the problem 
was started with the porous medium at constant temperature 
and oxygen concentration of 1100 degrees Fahrenheit and .006 

3 

lbm/ft , respectively. Figure V.7 shows the large cooling 

effect of the air entering the porous medium at x = 0 during 

the early part of the problem. As a result of this cooling, 

the heat generation began to raise the temperature interior 

to the porous medium before moving as a wave to the air inlet 

surface. The problem was terminated at 49 seconds with a 

porosity at x = 0 of .96. This was at the point of transition 

to the surface recession phase. Figure V.8 shows the behavior 

of the oxygen for this problem. At 49 seconds, the oxygen 

-4 3 

concentration at x = 0 was 4.1 * 10 lbm/ft and the penetra- 
tion of the oxygen was to x/L = .01. This again supports the 
assumption of treating the heat generation as a planar source 
and deleting the oxygen concentration equation during the 
surface recession phase. 

Other data generated by the model are spatial distributions 
of the following internal properties: particle diameter (d) , 

porosity (p) , specific permeability (m) , pressure (P) , pres- 
sure gradient, pore velocity (u) , Reynolds number (Re), heat 
transfer coefficient (h) , and surface area per unit volume (z) . 
As previously stated in the model development, these system 
properties are functions of temperature, and hence, change 
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FIGURE V.7 Sustained combustion produced by a high 

initial temperature. 
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with time. In addition, since surface recession may occur, 
the spatial thickness (L) of the porous medium is provided 
as output data. 

D. EFFECT OF GEOMETRIC PARAMETERS 

In the initial effort by Vatikiotis [2] , it was observed 
that, for a given set of boundary and initial conditions, the 
internal geometry of a porous medium significantly influenced 
combustion. Results were reported describing the dependence 
of combustion on particle size, and the thickness of the 
porous medium. As a refinement of that initial work, results 
are presented below which show the influence of the geometric 
parameters (i.e., permeability, porosity, porous medium thick- 
ness) on the minimum initial temperature needed to sustain 
combustion. "Minimum initial temperature" is defined here as 
the minimum constant-value temperature distribution (i.e., 
slope equal to zero) that will sustain combustion when used 
as an initial condition. A lower temperature will cause the 
porous medium to cool to ambient temperature. As will be 
shown in Section V.E, both the shape and magnitude of the 
temperature initial conditions influences the combustion proc- 
ess. In order to form a basis for comparison, the analysis 
was performed with uniform distributions for the initial con- 
ditions (the carbon and air having the same initial tempera- 
ture) . The initial condition for the oxygen concentration 
(also uniform) was taken as the concentration found in air at 
the initial condition temperature and ambient pressure. However, 
any oxygen concentration distribution, as an initial condition. 



would not have changed the results. The boundary conditions 
used in the analyses were the insulated boundaries, expressions 
III. 45 and III. 46, for the carbon temperature, and the Danck- 
werts conditions, expressions III. 47 through III. 50, for the 
air temperature and oxygen concentration. The results of the 
analysis were obtained by varying the initial conditions until 
the "minimum initial temperature" for sustaining combustion 
was obtained. A geometric parameter was then changed to a new 
value while all other parameters remained the same. A "mini- 
mum initial temperature" for this problem was obtained. This 
procedure was repeated for several values of the same geometric 
parameter while keeping all other parameters fixed. The re- 
sults obtained showed the "minimum initial temperature" which 
sustained combustion as a function of the respective parameter. 
The results for each geometric parameter are discussed 
individually . 

1 . Effects of Permeability 

Figure V.9 shows the dependence of the "minimum initial 
temperature" (defined previously in this section) on permea- 
bility. The description of the porous medium and the ambient 
conditions used in the following analysis is given in Table 
V.3. The range of permeability was determined by the magni- 
tude of the Reynolds number. The largest permeability 
— 8 2 

(m = .4x10 ft ) , at a pressure differential of 1.5 p.s.i., 
resulted in an average Reynolds number through the porous 
medium of approximately 5. It was felt that extending the 
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.2x10 .4x10 

FIGURE V . 9 — Permeability PERMEABILITY (ft 2 ) 

analysis. 



TABLE V.3 Geometry o-f porous medium and 

•for the permeability analysis 


ambient conditions 
in Section V.D. 



Particle Shape: 


spherical 


Particle Diameter (d >: 


(various) 


Unit Cell Thickness (D> : 


(various) 


Spatial Thickness o-f Porous Medium ([_)! 


1.0 in 


Porosi ty ( p ) 


. 476 


Permeability { m ) 


(various) 


Bulk Thermal Conductivity o-f Carbon ( k c ) ; 


86.0 Btu/ -f t-hr-F 


Bulk Specific Heat o-f Carbon ( Q ) : 


.231 Btu/lbm-F 


Bulk Density o-f Carbon <p c ) : 


70.3 lbrn/ft^ 


Thermal Emissivity o-f Particles (€ ): 


.9 


Ambient Temperature ( : 


80.0 deg-F 


Ambient Pressure ( P ) : 

CD 


14.7 psi 


Pressure Di -f -f erent i al (AP> : 


(various) 


Ambient Oxygen Concentrat i on (<b) : 


.0172 lbrn/ft^ 
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analysis to larger values of Reynolds number would exceed the 

limitations of the combustion model. Specifically, the 

assumption that Darcy's law governs the internal flow for 

Reynolds numbers greater than 5 would be questionable. Simi- 

-9 2 

larly, the smallest permeability (m = .162 x 10 ft ) , at a 
pressure differential of .1 p.s.i., gave an average Reynolds 
number of approximately .005. The lower limit Reynolds num- 
ber for which Darcy's law remains valid is not known. There 
is little information in the literature that addresses a lower 
limit Reynolds number. Although lower Reynolds numbers may 
be valid, the analysis was restricted to a minimum Reynolds 
number of .005. As can be seen in Figure V.9, the "minimum 
initial temperature" for sustaining combustion is a monotonically 
increasing function of decreasing slope of permeability. In 
other words, as permeability increases, the temperature needed 
for sustained combustion also increases. Important to the re- 
sults are that pore velocity also increases with increasing 
permeability. As stated previously, insulated boundaries 
were imposed on the porous solid. As will be shown in Section 
V.E, the results change when heat transfer occurs at the sur- 
faces of the porous solid. 

2 . Effects of Porosity 

Figure V.10 shows the dependence of the "minimum 
initial temperature" for sustaining combustion on porosity. 

The description of the porous medium and the ambient conditions 
used in the following analysis is shown in Table V.4. The 
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FIGURE V.10 --- Porosity POROSITY 

analysis . 



TABLE V.4 Geometry of porous medium and ambient conditions 

for the porosity analysis in Section V.D. 



Particle Shape: 


spherical 


Particle Diameter ( d ) : 


(various) 


Unit Cell Thickness <0 ): 


( var i ous) 


Spatial Thickness of Porous Medium <|_ >: 


1.0 in 


Porosi ty ( p ) 


(various) 


Permeability (m) 


(various) 


Bulk Thermal Conductivity of Carbon ( k c ) : 


86.0 Btu/f t-hr — F 


Bulk Specific Heat of Carbon (C c ): 


.231 Btu/lbm-F 


Bulk Density of Carbon (p c ): 


70.3 lbm/ft^ 


Thermal Emissivity of Particles (€ ): 


.9 


Ambient Temperature ( T^) : 


80.0 deg-F 


Ambient Pressure ( P ) : 

CO 


14.7 psi 


Pressure Differential (AP) • 


—.5 psi 


Ambient Oxygen Concentration : 


.0172 lbm/ft^ 
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results show that as porosity increases, the temperature 
needed for sustaining combustion decreases monotonically . 

This behavior is opposite, but less pronounced than that ob- 
served for the permeability. In addition, pore velocity de- 
creases as porosity increases. It seems reasonable that an 
increase in porosity provides more oxygen per volume of por- 
ous medium for combustion. The lowest value of porosity in 
the analysis (p = .48) was restricted by the idealized geometry 
shown in Figure III.l. For the idealized geometry, the low- 
est porosity occurs for the ratio of particle diameter to 
unit cell thickness, d/D, equal to 1. The maximum value of 
porosity in the analysis was limited by the highest porosity 
reported for spherical particles in a stable configuration. 

As stated by Scheidegger [9], this value is .875. 

3 . Effects of Porous Medium Thickness 

Figure V.ll shows the dependence of the "minimum 
initial temperature" for sustaining combustion on the porous 
medium thickness. The description of the porous medium and 
ambient conditions used in the analysis is shown in Table 
V.5. The range of thicknesses (.25"-4.0" ) investigated was 
limited (as in the case of permeability) by the Reynolds num- 
ber. As seen in Figure V.ll, the "minimum initial temperature" 
is a monotonically decreasing function with decreasing nega- 
tive slope of porous medium thickness. In other words, as 
the porous medium thickness decreases, the temperature needed 
to sustain combustion increases. In addition, the pore velocity 
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POROUS MEDIUM THICKNESS (in) 



TABLE V.5 Geometry of porous medium and 

■for the thickness analysis in 



Particle Shape: 

Particle Diameter (d ): 

Unit Cell Thickness (D): 

Spatial Thickness of Porous Medium ( L ): 
Porosity < p ) 

Permeability (m) 

Bulk Thermal Conductivity of Carbon ( k e ) : 
Bulk Specific Heat of Carbon ( c e ) : 

Bulk Density of Carbon ( p c ) : 

Thermal Emissivity of Particles <€ >: 
Ambient Temperature ( T^) : 

Ambient Pressure ( i^„) : 

Pressure Differential (Ap> s 
Ambient Oxygen Concentration ( <&) : 



ambient conditions 
Section V.D. 

spher i cal 
( various) 

( various) 
(various) 

. 476 

(various) 

86.0 Btu/f t-hr-F 
.231 Btu/lbm-F 

70.3 lbm/ft^ 

.9 

80.0 deg-F 
14.7 psi 

-.5 psi 
.0172 lbm/ft^ 
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decreases as the thickness increases. As in the case of the 
permeability, the results showed a wide range of temperatures. 
The "minimum initial temperature" for sustaining combustion 
may vary as much as 700 degrees Fahrenheit over a range of 
thicknesses for a particular porous medium. 

E. EFFECTS OF EXTERNAL PARAMETERS 

A similar analysis, as discussed in the previous section, 
was performed for the pressure differential. The boundary 
conditions for the carbon particle temperature equation were 
then changed from the insulated boundaries, expressions II I. 4 5 
and III. 46, to the radiation heat transfer boundary conditions 
expressions III. 51 and III. 52. This permitted heat transfer 
from the surfaces of the porous solid. The emissivity, e, was 
assumed equal to .9. With the new boundary conditions, the 
effects of permeability and pressure differential on the "mini 
mum initial temperature" for sustaining combustion were again 
analyzed . 

1. Effects of Pressure Differential 

Figure V.12 shows the effects of pressure differential 
on the "minimum initial temperature" for sustaining combustion 
The description of the porous medium and the ambient condi- 
tions used in the following analysis is shown in Table V.6. 

The range of pressure differentials considered (.1-2.0 p.s.i.) 
was governed by the Reynolds number for which Darcy's law re- 
mains applicable. This point was discussed in Section V.D. 

In the initial work of Vatikiotis [2] , the maximum pressure 
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FIGURE V.12 — Pressure PRESSURE DIFFERENTIAL (ps.i.) 

differential 
analysis . 



TABLE V.6 Geometry of porous medium and 

for the pressure differential 
Section V.E. 


ambient conditions 
analysis in 



Particle Shape: 


spher i cal 


Particle Diameter ( d ) : 


( vari ous) 


Unit Cell Thickness 


(various) 


Spatial Thickness of Porous Medium ( L )s 


1.0 in 


Porosi ty ( p ) 


.476 


Permeability ( m ) 


(various) 


Bulk Thermal Conductivity of Carbon ( k c > : 


86.0 Btu/f t-hr — F 


Bulk Specific Heat of Carbon ( C c ) : 


.231 Btu/lbm-F 


Bulk Density of Carbon ( £>, ) : 


70.3 lbm/ft^ 


Thermal Emissivity of Particles (€ >: 


.9 


Ambient Temperature ( T ) : 


80.0 deg-F 


Ambient Pressure ( P ) : 

CD 


14.7 psi 


Pressure Differential (AP) : 


(various) 


Ambient Oxygen Concentrati on (C jb^) : 


.0172 lbm/ft^ 
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differential considered was approximately .35 p.s.i. As seen 
in Figure V.12, the "minimum initial temperature" is a mono- 
tonically increasing function of pressure differential. The 
results show that, depending on pressure differential, the 
minimum temperature for sustaining combustion could vary as 
much as 400 degrees Fahrenheit for a particular porous medium. 
In addition, pore velocity increases as pressure differential 
increases . 

2 . Effects of Boundary Conditions 

In this section, the results of allowing heat transfer 
at the boundaries of the porous solid will be compared to the 
results obtained with insulated boundaries. Figures V.13 and 
V.14 show the results of the permeability and the pressure 
differential analyses, with and without heat transfer from 
the surface of the porous solid. The description of the por- 
ous medium and ambient conditions for the permeability and 
pressure differential analyses are shown in Tables V.3 and V.6, 
respectively. As stated previously, the radiation heat trans- 
fer boundary conditions were those shown by expressions III. 51 
and III. 52. As seen in Figures V.13 and V.14, heat transfer 
from the boundaries significantly affects the "minimum initial 
temperature" at the lower values of permeability and pressure 
differential (or at lower pore velocities) . The results show 
that when heat transfer occurs from the boundaries, the "mini- 
mum initial temperature" is no longer a mono tonic function of 
permeability and pressure differential. This behavior was 
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.5 1.0 1.5 2.0 

FIGURE V. 14 — Effects of PRESSURE DIFFERENTIAL (p.s.i.) 

boundary r 

conditions. 



suggested in the intial work of Varikiotis [2] , where in addi- 
tion to radiation heat transfer, convection heat transfer 
was considered at the boundaries. The additional convection 
heat transfer causes the effects produced by permitting radia- 
tion heat transfer at the boundaries to become more pronounced. 
Moreover, these results are consistent with those reported 
by Vulis [13] , Fontenot [51] , and Spalding [52] . The reason 
for this behavior was summarized in Section V.C and will be 
explained in Section V.G. 

3 . Effects of Initial Conditions 

The previous analyses of determining the effects of 
system parameters and boundary conditions were performed using 
constant-value initial conditions (i.e., slope of temperature 
distribution was equal to zero) . The constant-value tempera- 
ture defined the "minimum initial temperature" and was used 
as a basis for comparing the results. To also show the effects 
of the initial conditions on combustion, a single porous medium 
was subjected to numerous initial conditions. The analysis 
was performed as follows . The initial conditions of the car- 
bon and air temperatures were given a constant slope as shown 
in Figure V.15. With the slope fixed, the temperature distri- 
bution was shifted (i.e., moved up or down) until the thres- 
hold for sustained combustion was determined. That is, shift- 
ing the initial temperature distribution slightly lower would 
result in the porous medium cooling to ambient temperature. 

This procedure was repeated for several initial conditions of 
temperature having different constant slopes (both negative 
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FIGURE V.l 5 Constant-slope initial condition. 
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and positive) . The boundary conditions used in the analysis 
were the Danckwerts conditions, expressions III. 47 through 
III. 50, for the air temperature and oxygen concentration, and 
insulated boundaries, expressions III. 45 and III. 46, for the 
carbon temperature. The results of the analysis are shown in 
Figure V.16. In Figure V.16, the average temperature is plotted 
as a function of AT defined by T(L) - T(0) for the particular 
temperature distribution. The average temperature, T aV g* is 
defined as (T(L) + T(0))/2. As seen in Figure V.16, the curve 
is not symmetrical about AT = 0 . This result is atrributed to 
the direction of the pore velocity. The largest difference 
in temperature between the porous solid and the air (and great- 
est heat transfer) is at the air inlet region (i.e., x = 0) . 

At the air exit region (i.e., x = L) , the temperature differ- 
ence between the porous solid and the air is small. The com- 
bined effect of high heat transfer at the air inlet region 
and low heat transfer at the air exit region produces the 
following result. There is a greater likelihood that sustained 
combustion will occur for a positive slope temperature dis- 
tribution than for a negative slope temperature distribution 
having the same average temperature and AT as defined above. 

In other words, for a positive slope initial temperature dis- 
tribution, sustained combustion occurs at a higher temperature. 

F. EFFECTS OF PORE VELOCITY 

The way the parameters (i.e., permeability, porosity, 
pressure differential, and porous medium thickness) affected 
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TEMPERATURE DIFFERENCE (F) 



combustion is in proportion to their effect on pore velocity. 

In other words, the results of the previous sections showed 
that pore velocity was affected by the parameters in the same 
manner as the "minimum initial temperature" for sustained com- 
bustion (i.e., minimum temperature increased as pore velocity 
increased, and decreased as pore velocity decreased) . An 
exception to this occurred when heat transfer was permitted 
from the boundaries. This suggests that pore velocity is the 
dominant variable affecting the combustion process. Since 
pore velocity is a function of the parameters analyzed (i.e., 
permeability, porosity, pressure differential, porous medium 
thickness) , it seems reasonable that the "minimum initial 
temperature" for sustained combustion may be posed as a func- 
tion of pore velocity. With this approach, the "minimum initial 
temperature" for sustained combustion of a porous medium with 
unknown properties could be determined by specifying the 
magnitude of the pore velocity. Moreover, pore velocity can 
be represented by specific volumetric flowrate of filter velocity 
(i.e., V/A) which is easily measured. Figure V.17 shows the 
"minimum initial temperature", defined in Section V.D, as a 
function of the filter velocity at ambient temperature. The 
graph was made from the data obtained in the analyses of Sec- 
tions V.D and V.E (for insulated boundaries on the porous 
solid) . The data in Figure V.18 has a spread of approximately 
±100 degrees Fahrenheit. The equation for the curve shown 
in Figure V.17 is. 
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100 



(V.l) 



T t = 100 In 2 (Q/10) + 650 

where T gt represents the "minimum initial temperature" dis- 
tribution in degrees Fahrenheit, and Q is the filter velocity 
in ft/hr at an ambient temperature of 80 degrees Fahrenheit. 
Figure V.17 shows promise as a simple method for determining 
the "minimum initial temperature" for which sustained combus- 
tion can be expected in a porous medium. Perhaps by refining 
the formulation for calculating the pore velocity, there would 
be a smaller spread in the data. 

G. ANALYSIS OF RESULTS 

The Semenov model, discussed in Section III.C and illus- 
trated in Figure III. 2, may be used to explain the results 
reported in Sections V.B through V.F. Before continuing, cer- 
tain aspects of the Semenov model must be emphasized. First, 
the Semenov model does not govern combustion. It is an analy- 
tical tool which is used to describe the behavior of combustion. 
Secondly, the Semenov model as represented in Figure III. 2 
is for a point in the porous medium. A "Semenov surface" 
would better explain combustion in the porous medium. This 
will be illustrated below. It follows from the Semenov theory, 
as discussed by Vulis [13], that "ignition" conditions (i.e., 
sustained combustion) are determined by all the conditions of 
the combustion problem in a system. This is to say, that any 
change in system parameters (e.g., pressure differential, 
initial or boundary conditions) for a given porous medium 
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would result in a different Semenov surface. Lastly, the 
Semenov model is based on quasi-steady results. For example, 
the convection heat transfer, q^, and the heat generation, 

R^, will be equal at their intersection (see Figure III. 2). 
Conversely for the combustion model, at a particular point 

and time, q» may be smaller or larger than R depending on 

x, g 

whether the temperatures are increasing or decreasing. How- 
ever, if the rates of change in the temperatures are small, 
the Semenov model can be used to analyze the results of the 
combustion model. 

Figure V.18 illustrates a "Semenov surface" generated 
from the results of the last example in Section V.C (Figures 
V.7 and V.8). In Figure V.18, R is a heat generation sur- 
face, and represents the particle temperature-heat generation 
relationship during the problem. The q^ surface represents 
the convection heat transfer as a function of particle and 
air temperature (expression III. 14). The particular heat 
transfer surface in Figure V.18 represents the instant the 
point at x/L = 0 transitioned from the kinetic regime to 
the diffusion regime of combustion. This can be seen by the 
tangent point, I, between the heat generation surface and the 
heat transfer surface. A difference between "classical" 

Semenov S-curves (Figure III. 2) and those shown in the heat 
generation surface of Figure V.18 is the sudden drop in heat 
generation once a certain temperature has been reached. This 
is only observed in the region away from the air inlet (x/L = 0) 
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FIGURE V.18 Semenov surface. 
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and is attributed to the total depletion of oxygen in that 
region of the porous medium. One might say this is a weak- 
ness in the Semenov model when applied to porous media. 
Conversely, it might be argued that when the oxygen is totally 
depleted the combustion problem becomes a heat transfer prob- 
lem (in a pure sense) of which the Semenov model has no 
application. Also in the region where the oxygen has been 
depleted, it is not necessary for to intersect Rg. 

In order to visualize the results of the combustion model 
more clearly, the "Semenov surface" is abandoned for Semenov 
curves (illustrated in Figure III. 2). Figures V.19, V.20, 
and V.21 show that Semenov model representations for locations 
x/L = 0, x/L = .5, and x/L = 1, respectively, of example C in 
Section V.C (transient results shown in Figures V.7 and V.8). 

The heat generation curves are shown along with heat transfer 

% 

lines representing the film cooling at two different points 
in time. During this discussion, the term "film cooling" will 
be used to mean the heat transfer from the particle to the 
air. "Convection heat transfer" will refer to the transport 
of energy by the internal flow. In addition to film cooling, 
heat transfer by conduction and radiation occurs within the 
porous medium. This has an effect of slightly curving the 
heat transfer line upward (i.e., heat transfer becomes a non- 
linear vice a linear function of temperature) . However in 
this problem, calculations showed that the ratio of film 
cooling to conduction and radiation was approximately 5, and 
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therefore, conduction and radiation will not be considered 
in order to simplify the graphical analysis. 

Comparing Figures V.19-V.21 shows that the maximum heat 
generation rates are different depending on the location in 
the porous medium. It was mentioned in the general discussion 
of the Semenov model (Section III.C) that implicit in the 
shapes of the S-curves are the effects of oxygen concentra- 
tion. In other words, decreasing the oxygen supply will 
cause the S-curves to peak at lower temperatures. Since it is 
more difficult for the oxygen to penetrate as the depth of 
the porous medium increase (i.e., oxygen is consumed as the 
flow passes through the porous medium) , it is a reasonable 
outcome that the maximum heat generation rates decrease as 
the depth increases. The difficulty for the oxygen to pene- 
trate as the depth increases produces another effect as seen 
in Figure V.21. Unlike those in Figures V.19 and V.20, the 
S-curve and heat transfer curve does not have a tangent point 
which defines the "critical ignition condition". This be- 
havior is called noncritical combustion, and as discussed by 
Vulis [13] , is characteristic of oxygen poor combustion. 
Moreover, these observations are supported by the experimental 
results of Thomas, Stevenson, and Evans [15]. In their 
experiments, it was shown by decreasing the ambient partial 
pressure of oxygen, that a critical combustion process (i.e., 
one in which a "temperature jump" occurs) transforms to a 
noncritical combustion process. 
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As seen in Figure V.7, the temperature first rises in a 
local region of the porous medium (in this case near x/L = 1) 
before moving as a wave toward the air inlet surface. The 
direction of movement was a characteristic of all the prob- 
lems in which sustained combustion was observed. This behavior 
is explained as follows. Temperature rises in a region where 
the heat generation is greater than the heat transfer between 
the particles and the air. The region of higher temperature 
widens as heat is transferred by conduction and radiation. 

In addition, the increasing temperature simultaneously in- 
creases the reaction rate. At some moment, the oxygen enter- 
ing through the upstream side of the high temperature region 
is totally consumed within that region. When this occurs, 
the simultaneous increase in temperature and reaction rate 

will occur only on the upstream side of the region of high 

% 

temperature. As the increasing temperature moves towards the 
air inlet surface , the front of the region of depleted oxygen 
also moves nearer the air inlet surface. This results in the 
temperature response appearing to move in the form of a wave 
toward the air inlet surface. It is interesting to note, that 
the magnitude of heat transfer by conduction and radiation in 
the increasing temperature region was approximately twice 
that of the convection heat transfer at the outset of the 
example (Figures V.7 and V.8). It seems reasonable that by 
increasing the pore velocity to a sufficient magnitude, the 
region of increasing temperature would not propagate but blow 
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out (i.e., toward the x/L = 1 surface) . In fact, this is 
observed for the problems in which the porous medium cools 
to ambient temperature. 

Up to this point, the thermal response of the porous medium 
has been discussed in terms of "sustained combustion" or 
"cooling to ambient temperature" . The Semenov curves gener- 
ated from the results of the combustion model can be used to 
determine transition from kinetic to diffusion combustion. 

The kinetic and diffusion regimes of combustion are discussed 
in Section III.C. In a practical sense, the exothermic proc- 
ess occurring in the kinetic regime is synonymous with the 
term "oxidation". The process occurring in the diffusion 
regime is simply called "combustion". From everyday experi- 
ence, one knows that oxidation can be a slow process, occurring 
at low reaction rates. Conversely, combustion can produce 
considerable heat, reflected by high temperatures and reaction 
rates. As will be shown, transition from oxidation to combus- 
tion occurs rapidly. In addition, the porous medium may be 
undergoing oxidation in one part, and combustion in the other. 
It is important to understand the behavior of the thermal 
process at this level. Especially, if there is a need to 
closely control the process. 

Each point within the porous medium will transition from 
oxidation to combustion at different times. This can be 
seen in Figures V.19-V.21 which were obtained from the results 
of example C. The locations x/L = 1, x/L = .5, and x/L = 0 
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transitioned at 12.5 seconds, 22 seconds, and 41 seconds, 
respectively. Therefore at the start of the problem, the 
porous medium as a whole, was undergoing oxidation. More- 
over, transition from oxidation to combustion started at 
x/L = 1 and moved to x/L = 0 with time. This is consistent 
with the behavior observed for the temperature response. As 
discussed by Vulis [13] , the theoretical value for "critical 
ignition" temperature of the carbon can be calculated using 
the "N. N. Semenov equation" given by, 

T x = ^f-tl - (1 - 4R u T a /E) 1/2 ] (V . 2 ) 

where E is the activation energy of the Arrhenius expression, 

A 

R is the universal gas constant, and T is the absolute air 
u a 

temperature. This equation is derived by equating the 
Arrhenius expression III. 16 to the convection heat transfer 
expression III. 14, and requiring the slopes be equal. Note 
that the theoretical "critical ignition" temperature is a 
function of air temperature and activation energy only. Since 
the transition from oxidation to combustion at x/L = 1 shows 
a noncritical behavior, equation V.l does not apply. The 
theoretical "critical ignition" temperatures for locations 
x/L = 0 and x/L = .5 are 1831 degrees Fahrenheit and 1492 
degrees Fahrenheit, respectively. The values obtained by the 
combustion model (Figures V.19 and V.20) are approximately 
1825 degrees Fahrenheit for x/L = 0, and 1440 degrees Fahren- 
heit for x/L = .5. Once the "critical ignition" temperature 
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was reached, the increase to the "combustion" temperature 
occurred in less than 1 second for both cases . Examples A 
and B in Section V.C (Figures V.2-V.6) were the same problem as 
example C, discussed here, except combustion was initiated 
differently. Example C was initiated with the porous medium 
at a constant temperature known to sustain combustion. In 
examples A and B, the porous medium was subjected to a heat 
flux at the air inlet surface (x/L = 0) with the porous medium 
at ambient temperature. The amount of time for which the heat 
flux was applied determined if sustained combustion occurred. 
Sustained combustion was observed for example A in Section 
V.C. The amount of time the heat flux was applied in example 
B was not sufficient for sustaining combustion and the porous 
medium cooled to ambient temperature. It is interesting to 
note that the temperature at x/L = 0 in example A reached the 
"critical ignition" temperature of 1825 degrees Fahrenheit. 

The maximum temperature obtained at x/L = 0 in example B was 
approximately 1650 degrees Fahrenheit. There is good agreement 
between the theoretical values of "critical ignition" tempera- 
ture and those obtained by the combustion model . 

The pore velocity significantly affects the combustion 
process. This was discussed in Section V.F. Increasing pore 
velocity has two effects. First, the internal heat transfer 
coefficient, h, is increased, thereby, increasing the amount 
of film cooling. Secondly, the supply of oxygen is made greater 
(i.e., convection mass transfer of oxygen molecules increases). 
The combined effect is that by increasing pore velocity, higher 



112 



local temperatures are needed for transition of the porous 
solid from oxidation to combustion. In other words, assume 
that transition occurs at a particular air temperature for a 
given pore velocity. Increasing the pore velocity without an 
increase in air temperature may not be sufficient for transi- 
tion to occur. This is illustrated by Figure V.22. Increas- 
ing the pore velocity will increase the slope of the heat 
transfer line, and also change the shape of the heat genera- 
tion curve as shown. It is apparent that an increase in air 
temperature is needed to make the steeper-sloped heat transfer 
line tangent to the heat generation curve. In addition, higher 
pore velocities produce higher combustion temperatures once 
transition occurs. In the general discussion of the Semenov 
model, it was pointed out that the temperature dominates the 
combustion process in the kinetic regime. Increasing the 
oxygen supply will have the greatest effect in the diffusion 
regime (i.e. , increasing the combustion temperature) . This 
behavior described by the Semenov model, combined with the 
effects produced by convection heat transfer (energy trans- 
port by internal flow) and to a lesser extent, heat transfer 
by conduction and radiation, will determine if the porous 
medium will undergo sustained combustion for an increase in 
pore velocity. 

Figure V.23 shows the effects of decreasing the pore 
velocity. The slope of the heat transfer line becomes smaller, 
and the heat generation curve changes shape as shown. In 
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FIGURE V.22 Effects of increasing pore velocity. 
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FIGURE V.23 Effects of decreasing pore velocity. 
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particular, lower combustion temperatures will occur in the 
diffusion regime. Unlike the increasing pore velocity case 
above, there is less certainty in predicting the behavior of 
combustion when pore velocity is decreased. It appears that 
conduction and radiation heat transfer become more of a factor 
at the lower pore velocities. This seems reasonable since the 
ratio of conduction and radiation heat transfer to convection 
heat transfer increases with decreasing pore velocity. In 
addition, the effects of heat transfer from the boundaries at 
the lower pore velocities must also be considered (shown in 
Figure V.13 and V.14). The combined effects of all the heat 
transfer mechanisms, including the boundary conditions, will 
determine if the porous medium will undergo sustained combus- 
tion as the pore velocity is decreased. 

Lastly, the effects of the boundary conditions must be 
considered separately. The results of changing the boundary 
conditions were presented in Section V.E. It was observed 
that for heat transfer occurring at the surfaces of the porous 
solid (e.g., radiation boundary conditions), higher initial 
temperatures were needed to sustain combustion. Allowing heat 
transfer from the surfaces of the porous solid (either by 
convection or by radiation or by both) , affects the combustion 
process in the boundary condition regions as follows. Return- 
ing to the Semenov model, there is nothing explicitly asso- 
ciated with the boundary conditions that would change the shape 
of the heat generation curve. In addition, since pore velocity 
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determines the magnitude of the internal heat transfer coeffi- 
cient, it appears that the slope of the heat transfer line 
will also be unaffected. However up to this point, only film 
cooling was considered in the heat transfer line of the Semenov 
model. Noting that heat transfer from the porous solid boun- 
daries results in steeper temperature gradients in regions 
near the boundaries, heat transfer by conduction and radiation 
must also be considered. With the additional components of 
heat transfer, the heat transfer line will curve upward. This 
means that the heat transfer changes from a linear function 
to a nonlinear function of temperature. The change of the 
shape is illustrated in Figure V.24. As can be seen, the 
increased heat transfer by conduction and radiation may pre- 
vent transition from oxidation to combustion. This, in turn, 
will restrict the heat generation to the kinetic regime. It 
is reasonable to expect that as the reaction moves towards the 
air inlet boundary, the combined effects of greater heat trans- 
fer and lower heat generation will oppose sustained combustion. 
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HEAT FLUX 




FIGURE V.24 Effects of boundary conditions on 

heat transfer. 
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VI. CONCLUSIONS 



The main objective of developing a model which predicts 
the combustion behavior of a porous medium has been achieved. 
The results of the mathematical model are in good agreement 
with those obtained by experimental methods. Moreover, the 
analyses presented have shown that combustion and heat trans- 
fer in porous media is a complex process involving the inter- 
action of heat transfer and heat generation. This interaction 
depends upon the geometric parameters of a porous medium, 
boundary conditions, environmental and initial conditions. 

The behavior of a system can change radically by altering any 
number of parameters. This point was demonstrated by examples 
A and B in Section V.C where a difference of one second in 
applying a surface heat flux meant combustion or extinguish- 
ment. The results have also shown that conduction and radia- 
tion between particles may play a significant role in deter- 
mining combustion behavior and should be accounted for. The 
computer program makes it possible to look at a large number 
of cases. Similar analyses by experimental methods would be 
economically impractical. 

Based on literature surveys performed during this inves- 
tigation, it appears that much of the engineering develop- 
ment associated with porous media involves a trial and error 
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process. It seems reasonable that an analytical model would 
be an essential tool for the engineer to employ in the de- 
sign process. In this investigation, the model has been 
assessed for its applicability, and to some extent, its 
accuracy. It is hoped that in the future the model will be 
used to determine the effects of the design variables (i.e. , 
the geometric parameters, pressure differential, etc.) on 
performance. Specifically, the model shows promise for 
evaluating the combustion efficiency and stability of a sys- 
tem. It is in this capacity that the combustion and heat 
transfer model will be of greatest value. 
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APPENDIX A 



FORMULATION OF FIELD EQUATIONS 



1. PRESSURE DISTRIBUTION EQUATION 



Darcy's law for one dimensional flow, neglecting body 
forces, is, 



Q 



m ,dP> 
y dx 



( A. 1 ) 



Substituting in the Dupuit-Forcheimer relation, and solving 
for u, equation A.l becomes. 



u 



m ,dP, 
py 'dx 



(A. 2) 



The continuity equation (derived in Appendix B) is. 



3 (PP a ) 
3t 




u) 



0 



(A. 3) 



Substituting equation A. 2 into equation A. 3 yields. 



3 _ 3_.^a dP. 
3t 3x' y dx 



(A. 4) 



Expanding terms, equation A. 4 becomes, 



d 2 P _1_ ^a + l3m_l3ysdP y_ 8(pp a } 

^ x 2 'p a 3x m 3x y 3x dx mp^ 3t 
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Equation A. 5 with associated boundary conditions (presented in 
Section II. B) is solved for the pressure and pressure gradient 
distribution. The pressure gradient is then substituted into 
equation A. 2 to obtain the pore velocity. 

2. POROUS SOLID HEAT TRANSFER EQUATION 

To perform energy balances on both the porous solid and 
on the air, a differential volume of porous medium was segre- 
gated into respective volumes of constituents, that is, 
dV = (l-p)dV for the solid, and dV = pdV for the air (shown 
in Figure A.l). The convention used for the energy balance 
of an arbitrary differential volume, dV, is, 

Heat into Heat _ Heat out + Increase in 

dV Generation of dV Internal Energy 

The heat transfer mechanisms considered for the carbon 
particles are conduction, radiation heat transfer between 
particles, convection heat transfer from the particles to 
the air, and heat generation. Applying the above convention, 
the energy balance on a differential volume of porous solid 



is 



(1 -P | '3 C ond dA 



X 




(A. 6) 



(1 -P ), 3cond dA | . + (l-P)5 ra d dA | _ + q 

x+dx 1 x+dx ' 



conv 



dA' 
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POROUS MEDIUM AIR SOLID 




FIGURE A.l Separating a differential volume of 

porous medium into respective volumes 
of solid and air. 
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Representing terms on the right side of expression A. 6 by 
Taylor series expansions (neglecting higher order terms) , 
the energy balance becomes , 



< 1 -e ) ‘J c ond dA | + ( 1 ' p)g rad dA | + Sen 3 *' 

X X 



(A. 7) 



“ ( 1 -P»lcond dA | + u‘ ( 1 'P> g cond ldxdA + ( 1 -P> g rad dA | 

X 



X 



+ 3x [ (1 “P )c 3>. a ^] dxdA + <3™^’ + (1-P)q^<3v 



'rad 



conv 



l int 



Subtracting terms, and rearranging, expression A. 7 becomes. 



- k l( 1 -P )g cond )dV - fe l( 1 -P )g rad IdV ‘ '’conv^' 



+ q dA ' = (l-p)q. , dV 

^gen ^ tnt 



(A. 8) 



Substituting the following expressions into equation A. 8, 



3T 



'cond 



= - k 



e 3x 



Fourier's law 



(A. 9) 



3T 



'rad 



= - k 



r 3x 



Radiation analogy to (A. 10) 

Fourier's law 



q conv = h(T c ” T a ) Newton's law 



(A. 11) 



q = R 

^gen g 



Heat generation 



(A. 12) 
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Internal energy 



(A. 13) 




yields 




] dV - h (T -T ) dA ' + R dA 1 
c a g 



(A. 14) 



(1_p)p c c c — dV 




Dividing through by dV, and defining cLA'/ciV as z, the specific 
internal area (i.e., surface area per unit volume) , equation 
A. 14 becomes. 



The expressions used to obtain the values of the properties 
and parameters in equation A. 15 are presented in Section 
III .E. 

3. AIR HEAT TRANSFER EQUATION 

The formulation of the air heat transfer equation will 
begin with the general one dimensional energy equation. 



I— [ (1-p) (k +k ) 
9x * e r 



3T 

] - hz(T -T ) + R z 



(A. 15) 





3 _ 

3x 



3T 

-r-^-) + hz (T -T ) + UF 



(A. 16) 
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As in the porous solid heat transfer equation, the time and 
position dependent porosity appears inside the differential. 
Expanding terms and neglecting body forces, the air energy 
equation becomes. 



De 

P p aDt + P p a 



D(|u 2 ) 

Dt 



3 T 

+ hz( V T a> 



(A. 17) 



- fx (p u P * ' k (puT xx> 



Consider the momentum equation for the x-direction (neglecting 
body forces) , 



3t (pp a u) = 



3 , 2 , 3 , . 3 , 

' a5 <pp a u > - pT xx - ^x pP 



(A. 18) 



and the continuity equation. 



!t‘ pp a> “ 



- u fc (pp a> 



3u 

p p a 3x 



(A. 19) 



Multiplying the continuity equation through by u, and substi- 
tuting this into equation A. 18, the momentum equation becomes, 



3u 

p p a 3t 



3 U 3 / \ 3 / X 

- pp a u H - 5S (pP) - H (pT xx> 



(A. 20) 



Multiplying equation A. 20 by u and noting that. 



D U 

Upp a Dt 



3u , 2 3u 

upp a 3t + pp a u 3S 



(A. 21) 
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equation A. 20 becomes. 



u p 



A Du 
p a Dt 



- “ fx<P P > - U k ( P T XK> 



3x 



(A. 22) 



The energy equation A. 17, after substituting 
1 2 

DC-^-u ) /Dt = uDu/Dt and expanding terms, is. 



De 



Du 



P°aDt + P p a u Dt 



31 T 

= + hz (T-T ) 

dx a dx c a 



P P - a iiESl- 
r 3x 3x 



3 , . 3u 

U 3x P T xx } P T xx 3x 



(A. 23) 



Substituting equation A. 22 into the above energy equation and 
cancelling terms, equation A. 23 becomes. 



De 

P p a Dt = 



3^ (pk a 



3T 

a 

3x 



+ hz (T -T ) - pP^ - 
c a ^ 3x 



3u 

^ xx 3 x 



(A. 24) 



The viscous dissipation term, px 3u/3x, is neglected because 

XX 

the fluid is a gas flowing at a low velocity. Therefore, the 
energy equation for the air in the porous medium is. 



De 

p p a Dt 



a 3T a 

_ (pkaT ^) + h Z (T c -T a ) - PP^ 



(A. 25) 



With specific enthalpy for a gas defined by. 



k = e + P/p 



(A. 26) 
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pe/pt can be expressed as. 



pe _ Dk 1 DP 
pt pt ” p a pt 

a 



Multiplying the continuity 



P Dp a 
2 Pt 



(A. 27) 



equation through by P/(pp ) gives, 

3 . 



P Dp a = _ P PP P 9u 

7 ? Dt “ pp Pt “ p 3x 

P a a 



(A. 28) 



Substitution of the above expression into expression A. 27, 
the substantial derivative of internal energy can be expressed 
as. 



pe _ pfr_ _1_ pP _ _P_ 3_u P PP 

pt Pt”p_pt”p 3x pp pt 

da. a. 



(A. 29) 



Substituting expression A. 29 into equation A. 25, the energy 
equation reduces to. 



P P 



D h. 

a pt 






ST 

— -) 
a 3x ' 



+ hz(T -T ) 
c a 



(A. 30) 



The specific enthalpy can be represented in terms of tempera- 
ture , as , 



dh = T ds + - dP (A. 31) 

P 

where s is specific entropy. For a perfect gas, ds may be 
expressed by. 



128 



ds 



(A. 32) 



■ C P ¥ - £ dp 



Substituting expression A. 32 into expression A. 31, and can- 
celling terms, enthalpy can now be written as. 



dh = c dT 
P 



(A. 33) 



or , 



Dh 

Dt 



dt_ 



= c 



a Dt 



(A. 34) 



Substituting expression A. 34 into equation A. 30, the air 
energy equation or heat transfer equation becomes, 



a 3T 

3* <p k a TT> 



3T 

cl 

p 0 U C -r 

^ ^a a 3x 



+ hz (T -T ) 
c a 



(A. 35) 



+ 5t ( P P) - 



P p a c a 



3T 

a 

at 



Initial results showed that pP changed slowly with time. 
Therefore, the substantial derivative of pP, as shown here. 



Bt ( P p> = ft (pP) + u fe (pP) 



(A. 36) 



is reduced to u3(pP)/9x. The expressions used to obtain the 
properties and parameters in the coefficients of equation 
A. 35 are presented in Section III.E. 
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4. OXYGEN MOLECULE DIFFUSION EQUATION 

The final consideration in formulating the field equa- 
tions for the model is the transport of oxygen molecules. 

The oxygen molecule transport equation is obtained by a con- 
servation of species belance on the differential volume of 
air, dV = pdV. The convention used for the species balance 

cL 

is given by. 



C >2 into _ O 2 out + O 2 + O 2 

dV of dV Consumption Accumulation 

The transport mechanisms considered were diffusion resulting 
from concentration gradients, convection, and consumption of 
oxygen by combustion. Applying the above convention, the 
species balance on the oxygen becomes. 



P m. • dA + pm dA = pm, . -,dA 
r diff ^ conv r diff 1 

X 1 X 



(A. 37) 



x+dx 



+ p m dA 
v conv 



x+dx 



+ m dA ' + p m dV 

cons * acc 



Representing the terms on the right side by Taylor series 
expansions (neglecting higher order terms) , the species balance 
becomes , 



P” di ff dA | + P m conv dA | * P ra diff dA 

X 1 X 



(A. 38) 



x 



+ tt— ( pm, . -j.) dxdA + pm dA + ^— (pm ) dxdA + m dA' +pm dV 
3x'^ diff ^ conv 3x * conv cons ^ acc 

x 



130 



Cancelling terms and rearranging, equation A. 38 becomes. 



^ — (pm , . dV - ^ — (pm ) dV - m dA' 

3x ^ diff 3x y conv cons 



pm acc dV (A<39) 



Substituting the following expressions into equation A. 39, 



m diff 



_ p 14. 

e 3x 



Fick ' s law 



(A. 40) 



m = u d> 

conv y 



Convection transport (A. 41) 



m. 



cons 



= R, 



Consumption by 
combustion 



(A. 42) 



m 



acc 



= M. 

3 1 



Accumulation 



(A. 43) 



yields. 



|j(p P e |i)dV - |j(up*)dv - R^dA' 



P It dV 



(A. 44) 



Dividing both sides by dV, and letting dA'/dV equal the 
specific internal area, z, the oxygen molecule diffusion 
equation becomes , 



H ( e p eH ) - H (u p4’> - r o 2 z 



pH 



(A. 45) 



The methods and expressions for obtaining the properties and 
parameters in the coefficients of equation A. 45 are presented 
in Section II. F. 
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5. TRANSFORMATION OF FIELD EQUATIONS FROM A FIXED COORDINATE 

TO A MOVING COORDINATE SYSTEM 

As discussed previously, the reaction rate at the air in- 
let surface of the porous medium results in surface recession. 
To account for this surface recession, the field equations 
must be transformed from a fixed coordinate to a moving coor- 
dinate system. The method of transforming the field equations 
will be shown for only the porous solid heat transfer equa- 
tion since the approach for the other field equations (i.e., 
air heat transfer, and combined Darcy's law and continuity 
equation) is identical. 

The porous solid heat transfer equation with the modified 
heat generation term is. 



9 T 

btd-pxw-j # 1 - hz( W + v A(x = 0) 



3T 



d*P)P c c cTt: 



(A. 46) 



Since the x coordinate is a function of time during surface 
recession (e.g., T (x(t),t)), the time derivative term in 
equation A. 4 6 must be expanded using the chain rule. 



fp[T (x(t) ,t) ] 

d t C 



9T 

= ( — ) 

v 3x ' 



<at>. 

t l 



3T 

+ <TT> 



X 



(A. 47) 



The x coordinate in the field equations is nondimensionalized 
by n = x/L. Substituting n for x, expression A. 47 becomes. 
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X 



(A. 48) 



3 T 

( — — ) 
v 3x ; 



(— ) 
t dt i 



3 T 

+ <TmT> - < 






3 T 

+ <TF> 



Noting that. 



3jl _ 1 

3x L 



(A. 49) 



and. 



dx _ dn dL 

dt - L dt + n dt 



= n 



dL 

dt 



expression A. 48 becomes, 



3T 

( — — ) 
v 3x ' 



(— ) 

. l dt J . 
t l 



3T 

+ <TF J 



j t 3T 

H ( dL > c 

L v dt ; 3n 



3T 

c 
3 1 



(A. 50) 



(A. 51) 



Upon substituting n = x/L into the remainder of equation A. 46, 
and substituting expression A. 51 for the time derivative term, 
the porous solid heat transfer equation for the surface 
recession problem is. 



3 T 

L " 2 |-[(1-P) (k +k )— f] - hz ( T -T ) + R • 6 A (r, = 0) 
3r) e r 3n c a g 



. 3T 

= (1-p) p c [2-L-s-S. + 

c c L 3n 



3T 

c 

3t 



(A. 52) 



As stated in Section III.G, the thickness, L, as a function 
of time, and L can be obtained from expression III. 44. The 
air temperature and combined Darcy's law and continuity 
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equation appear similar to equation A. 52 for the surface 
recession problem. 
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APPENDIX B 



FORMULATION OF AUXILLIARY EQUATIONS 

1. CONTINUITY EQUATION 

The law of conservation of mass requires that the substan- 
tial derivative of the fluid mass in a differential volume be 
zero. Therefore, the continuity equation for a fluid in a 
porous medium is expressed by, 

^<PP a dV) = 0 (B.l) 

or in an equivalent form, 

b (E>p a> + u fe ( P p a> + P p als = 0 tB ‘ 2) 

2. RADIATION HEAT TRANSFER ANALOGY TO FOURIER'S LAW 
Radiation heat transfer in the model was represented by 

an analog to Fourier's law of conduction heat transfer shown 
by. 



Orad = - k r Tx (B ' 3) 

where is an equivalent conductivity of the particles due 
to radiation. The development is as follows. Assuming air 
to be transparent to radiation, and treating the idealized 
geometry of the porous medium (Figure III.l) as a series of 
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closely spaced walls, the net radiation heat flux between 
adjacent walls is. 



*rad 



- — — (T^ - ) 

2-e^x x+dx' 



(B . 4 ) 



where and T X+ £ X are the respective absolute wall tempera- 
tures, e is the emissivity of the carbon, and a is the Stefan- 
Boltzman constant. Expanding <3 ra £ in a Taylor series about 

A 

T x+dx' an< ^ ne glecting higher order terms, the series expan- 
sion may be written as, 



q = — — — (T^ - T^) - ^ a£ dx 

q rad 2-e ' i x x ; 2-e x+dx 3x ax 



(B . 5) 



Simplifying, the above expression becomes. 



*rad 



4ae p 



3T 



2-e x+dx 3x 



dx 



(B. 6) 



Equating expressions B.3 and B.6, and noting that 9T/3x = 3T/3x, 

k becomes, 
r 



k 



r 



4ae 

2-e 



dx T 



x+dx 



(B . 7) 



where dx is now equal to 5, the pore diameter. From the close 
spacing of the particles, the temperature difference will be 
small as compared to the magnitude of the temperature. Noting 
this, the average absolute temperature of the carbon particle 
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The equivalent radiation con- 



may be substituted for T x+C j x * 
ductivity expression becomes. 



k 



r 



4ae<S ~3 
2-e c 



(B . 8 ) 



and the radiation heat transfer from particle to particle 
may be represented by. 



*^rad 



4cre<5 3,3 ^ T c 
2-e c 3 x 



(B . 9 ) 



For small pore diameters equation B.9 will be a good approxi- 
mation of the radiation heat transfer between particles. 



3. POLYNOMIAL APPROXIMATIONS OF THERMAL PROPERTIES 

Relations giving the dynamic viscosity, thermal conduc- 
tivity, and specific heat at constant pressure of air at 
different temperatures were required. A simple method to 
obtain values for these properties is to fit empirical data 
with 2nd order Lagrange polynomials. 

The general form of the 2nd order Lagrange plynomial is. 



(T i -T 2 ) (T i -T 3 ) (T i" T 3 ) ( T i -'V 

k i = <VT 2 > < Tl — t 3 ) k i + <T 2 -T ) <T 2 -T ) 



(T.-T, ) (T . -T„) 

_i_ 1 1 ^ ]r 

(T,-T, ) (T,-TJ 3 



(B.10) 



where k^ is the property value at the ith temperature, T^ . 
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Choosing three temperatures that are representative of 
those observed during the analysis, the corresponding values 
of the properties are. 



(subscript) 


Temp . 


y 


k 

a 


c 

a 




deg-F 


lbm/f t-hr 


Btu/f t-hr-F 


Btu/lbm-F 


1 


80 


4.7973xl0 -2 


1.5161xl0 -2 


.24020 


2 


1340 


1.0045xl0 _1 


3.9013xl0 -2 


.27268 


3 


3140 


1.5725xl0 _1 


7.1647xl0 -2 


.31196 



Applying expression B.10 to each set of properties results 
in the following set of polynomials, 

u = -3.308 X 10 “ 9 T 2 + 4.633 x10~ 5 T +4.427xlo" 2 (B.ll) 

cl cl 

k = -2.608 x 10‘ 10 T 2 + 1.930 x 10" 5 T + 1.361xl0 _2 (B.12) 

a a a 

c = -1.29 3 x 10~ 9 T 2 + 2.758 x 10 -5 T + .238 (B.13) 

a a a 

Each expression gives property values within two percent of 
the data for temperatures to 3000 degrees Fahrenheit. 

4. JUSTIFICATION OF THE DANCKWERT ' S BOUNDARY CONDITIONS 
The Danckwerts ' boundary conditions for the air heat 
transfer are. 



3T 

a 

3x 



= P. 



u c (T - T ) 
a a °° 



( B . 1 4 ) 
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0 



(B . 15 ) 




Bischoff [41] presents a discussion of these boundary condi- 
tions applied to mass diffusion equations in porous media. 

An analogous discussion is made here for the Danckwerts' 
boundary conditions applied to the air heat transfer equation. 

The porous medium with entrance and exit regions is shown 
in Figure B.l. The air temperature, T , for each section of 

cl 

the region is distinguished by subscripts, as are the proper- 
ties. The air heat transfer equations for each of the regions 
are as follows. 



d 2 T 



dT 



dx 



- p u, c 
M a 1 a 



dx 



= 0 



x < 0 



(B . 16) 



- dT dT 

dx (p k a =asf> ' °a uc a dF 



+ hz (T -T ) = 0 0 < x < L (B . 17 ) 

C 3. — — 



d 2 T 



dT 



‘a , 2 
dx 



- p u, c 



a 2 a dx 



= 0 



x > L 



(B . 18 ) 



with the boundary conditions. 



T (-oo) = T (B . 19 ) 

a l 

T (0) = T (0) (B . 20 ) 

a. a 
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Flow e£> 




FIGURE B.l Justification of the Dankwerts 

boundary conditions. 



140 



(B . 21) 



p a u l c a T a 1 (0) ' k alf [T a 1 (0)! = 



pp a uc a T a <0) ' pk aS [T a <0)! 



T (L) = T <L) 

d n d 



(B . 22) 



p u~c T 
K a 2 a a. 



(1 > - k a!*‘V L) > ‘ pp a uc a T a (0 > - pk a|-tT (L) ] 



T (°o) = finite 

a 2 



(B . 23 ) 
(B . 24) 



The above boundary conditions impose the restriction that 
there is no convection heat transfer from the porous solid 
to the air at the boundary surfaces. The difficulty of mixing 
Danckwerts' and convection heat transfer boundary conditions 
is discussed in Section III.H. For convenience, the thermo- 
physical properties will be treated as constant in the entrance 
and exit regions. 

An analytical closed-form solution of this set of equations 
is unknown because of the nonlinearity of the temperature 
dependent properties in equation B.17. However, the solutions 
of equations B.16 and B.18 are. 



T 

a 



1 



K 



+ K 2 exp ( 



p a U 1 C a 



x) 



x < 0 (B . 25) 



T 

a 



2 



P a U 2 C a 

k 3 + k 4 ex P (- ^k ” x) 

a 



x > L (B . 26) 
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Applying boundary conditions B.19 and B.24, the following 
results are obtained, 



K 



1 



T 



CO 



(B. 27 ) 



K 



4 



= 0 



(B. 28) 



The solutions, B.25 and B.26, become. 




T + K- 

oo 2 



p U, C 

exp — — x) 



x < 0 (B . 29 ) 



T = K_ x > L (B . 30) 

It would be necessary to have the solution for the nonlinear 
equation B.17 to solve for and . However, the constants 

need not be known to continue with the analysis. 

From equation B.29, 



T (0) = T + K 9 

a, «> 2 



(B . 31) 



and , 



d_ 

dx 




(0) ] 



p u. c 
H a 1 a 



K, 



Substituting these into equation B.21 gives. 



(B . 32) 
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p a u l c a T °° + p a u l c a K 2 - p a u l c a K 2 “ P p a uc a T a (0) -P k a3* [T a (0)1 

(B . 33) 

Cancelling terms and noting = up (i.e., u^p^ = u 2 P 2 ) 
yields , 



p p u c T 
^ M a a 00 



- P p a uc a T a (0 > ' P k ad5>V 0)1 (B ' 34) 



Rearranging, the Danckwerts 1 boundary condition at x = 0 is. 



k a $? [T a (0)1 



p a uc a (T a ' T »> 



(B . 35) 



Noting that is a constant, the derivative of equation B.30 
is. 



4 - T 

dx a. 



= 0 



(B.36) 



and. 



dJ tT a 2 (L)1 = 0 



(B. 37) 



Sbustituting these expressions and equation B.22 into equa- 
tion B.23, and noting u 2 = up (i.e., u^p^ = u 2 p 2 ) 9 ^- ves ' 



P p a uc a T a (L) - k a l? [T a (L)1 = 



P p a uc a T a (L) 



(B . 38) 



Upon cancelling terms, the second Danckwerts' boundary 
condition becomes. 
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0 



(B . 39 ) 



d_ 

dx 



[T, (L)] = 

d 



An important consideration for using the Danckwerts ' 
boundary conditions is that it simplifies the analysis since 
equations similar to B.17 may be solved independently without 
having to consider entrance and exit regions. 

5. PRESSURE DIFFERENTIAL CALCULATION RESULTING FROM EXTERNAL 
FLOW AT SURFACE X = L 

To obtain the pressure differential, AP, for simulating 
the conditions of Fontenot's experiments [1] (discussed in 
Section III.H) , Bernoulli's equation was used. The following 
observations showed this to be a valid assumption. Schlichting 
[54] points out that for the ratio of u/U^ in the range of 
.0001 to .01, the effects of "blowing" or "suction" on the 
potential flow over the external surface of the porous medium 
may be neglected. A typical value of u/U^ for the model at 

which U = 25 knots was .0028. For steady flow over a flat 

00 

plate, the flow field outside the boundary layer may be 
described by Bernoulli's equation. This is a direct result 
of the Navier-Stokes equation. In addition, the pressure 
gradient across the boundary layer may be taken as zero. 
Therefore, Bernoulli's equation, 

2 

f — + = constant (B.40) 

’ p a 2 

may be used to obtain the pressure differential across the 
plate. The parameters in equation B.40 are defined as follows: 
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P, pressure; p , density of air; U, free stream velocity of 
the air over an external surface of the porous medium. For 
the model, the density of the air was approximated by the 
ideal gas law as. 



= P/R. 



a 



(B . 41) 



where R is the gas constant of air, T is the absolute tern- 
ci a. 

perature of the air, and P is the pressure. Substituting 
equation B.41 into equation B.40 gives, 



/ 



dP + 




constant 



(B . 42) 



Upon integrating, equation B.42 becomes, 



R a T In (P) 
a a 




constant 



(B . 43) 



or, 



R 



a 




u; 



ln(P ] _) 




ln(P 2 ) + 




(B . 44) 



Letting , 




and substituting these into equation B.44, yields, 
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/V 



(B . 45 ) 



R 

a 



T 

00 



ln(P) 

00 



R T In (P ) 
a 00 L 




Rearranging, equation B.45 becomes. 



R 



a 



/\ 




ln (Pco/P L ) 




(B . 46 ) 



Taking the exponential of both sides, and solving for P , 

Li 

results in. 



= exp (-U“/2R 



T ) 

OO 



(B . 47) 



From the above expression, and noting that AP = P -P , AP 

Li 00 

may be expressed as. 



AP = P [exp(-uf/2 R T ) - 1] 

oo oo 3. 00 



(B. 48) 



Expression B.48 was used to approximate the pressure differ- 
ential across the porous medium when simulating the conditions 
of Fontenot's experiments [1], 
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APPENDIX C 



NUMERICAL FORMULATION 
1. FINITE ELEMENT METHOD 

The solution of the system of nonlinear, coupled, partial 
differential equations given by equations III. 32, III. 38, and 
III. 39, subject to boundary and initial conditions, was ob- 
tained by a Galerkin formulation of the finite element method. 
The solution of equation III. 12 was obtained using a shooting 
method . 

a. Galerkin Formulation 

A Galerkin formulation of the Finite Element Method 
was used to obtain solutions of the porous solid and air 
energy equations, and the oxygen diffusion equation. A con- 
venient form of equations III. 32, III. 38, and III. 39 was 
used in the formulation as shown by. 



3 T 

l" 2 |— [ (l-p) (k +k )-r^] - hz (T -T ) + R z 
9ri e r 3 n c a g 



3T 



= (l-p) pc 



c c at 

(C.l) 



— 2 3 -1 

L i_(p k _f) - p p c uL 1 hz(T -T ) 

dri a 3ri sl 3 l 3f) c a 



(C. 2) 



+ uirl !^e p) 



3T 



= PP c 



a a 3t 



L ' 2 k'P V e fl 1 - L_1 frT <p u ♦> - R 0 2 - P II 



(C . 3 ) 
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where the spatial coordinate, x, was nondimensionalized by 
n = x/L. 

The closed domain (0,1) was partitioned into (n-1) 
contiguous elements of variable length £^, i = l,...,n-l. 
This defines an n nodal point model. The three field varia- 
bles, T , T , 4> were approximated by, 

C cl 



T c (n,t) 


= ^(rut) 


= l G(n)0 1 (t) 


(C . 4 ) 


T a (n,t) 


= 4^2 (H/t) 


= l G(n)0 2 (t) 


(C . 5) 


<Mn #t) 


= 4< 3 (n,t) 


l G(n)0 3 (t) 


(C . 6 ) 


where G^ , for i = l,...,n 
tions with local support. 


is a set of specified 
and the sets { 9 ^, 92 / 0 


basis func- 

3 , i 1 , • • 



are the solution coefficients to be determined. The G. were 

1 



selected to satisfy the condition G^(rij) = 6 ^ ^ , where the 

Jc Ic 

Kronecker delta, 6 . ., is defined by 5 . . = 1 for i = j, and 

1 J 1 j 

5^^ = 0 for i 7 ^ j . As a result, 9^, a nd 9 3 are the 

values 4 ^' ^3 at nodal points (i.e., 9^ (t) = 4 ) j_ (n^ • 1) ) . 

1 

Linear interpolation functions (shown in Figure C.l) 
were used as the basis functions. These are the lowest 
polynomial functions which provide the necessary function 
continuity. 

As a measure of error, a residual function, r^, is 
defined for each field equation by, 
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FIGURE C.l Linear shape functions used in the 

Galerkin formulation of the FEM. 
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1 



1,2,3 



(C . 7 ) 



r i (n , t) = A i (ip ) “ 

where denotes the spatial operator of the ith equation. 

For convenience, the following convention for differentiation 
is adopted. 




(C. 8) 



2-y ( ) = ( )" (C . 9 ) 

3n 




(C.10) 



For field equations C.l, C.2, and C.3, the residuals are. 



- n n 

r. = lf 2 [(l-p) (k +k ) l G J 0 . ]'-hz l G. (0, -0, ) (C.ll) 

1 r i=l 1 1 1 i=l 1 i i 



+ R 0 

g 




n 



I G i9 

i=l 1 



<1-P)P c = c 



n 

l G.0, 
iii 1 



-9 n -1 n 

r 2 = L (pk l G 0 - )' - pp c uL 1 l G!0 2 (C . 12 ) 

2 a i=l 1 2 i a a i=l 1 2 i 



n n 

+ hz l G. (0 -0, ) + uL i (pP) '-PP c l G. 0 2 

i=l 1 X i 2 i a a i=l 1 2 i 
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(C . 13) 



• - L" 2 (pP e l G!0 )' - L _1 (p U l G.0 )' 

i=l i i=l i 



- R 0 " 1 z l G.0. -p T G.0, 

2 i i=l l i=l l 



n 



n 



where the coefficients multiplying the response variables are 
themselves functions of the response variables, and thus, 
the equations are nonlinear. In accordance with the Galerkin 
method, the final system of ordinary differential equations 
was obtained by setting each residual, r ^ , orthogonal to each 
basis function, G^ , that is, 

1 l _ 1 ,2 , . . . , n 

/ G. r . dn = 0 (C . 14 ) 

0 1 J j = 1,2,3 

The 3n ordinary differential equations given by equations C.14 
retain the character of the original set of partial differ- 
ential equations. Thus, linear field operators transform to 
matrix operators and nonlinear, coupled field operators become 
nonlinear, coupled algebraic operators. Incorporation of 
the boundary conditions resulted in 3n nonlinear coupled 
ordinary differential equations, 

A ( t) 4;(0 1 ,0 2 ,0 3 ) + F ( t) = B(t) ifj (C.15) 

subject to initial conditions, where B is a 3n >: 3n matrix. 
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u ;> 



is the operator associated with the field operator in 
expression C.7, and F is an excitation vector. 

Adopting the convention. 



<G l > { 9 i > 



n 

l 

i=l 



G i 9 i 



(C . 16 ) 



and applying the operation of expression C.14 with an integra- 
tion by parts on the second order derivatives gives. 



(G i }L" 2 (l-p) (k e +k r ) <G ± > ' (9 I J 



L -2 ( 1-p) (k e +k r ) 



1 

/ { G . } ' <G . > ' dri 
0 3 



1 



-hz / 
0 



{G i }<G j >dn(9 1 } 



+ hz 



1 

/ {G. } <G . >dri { 9 _ } 
0 1 J *■ 



(C . 17) 



+ R z9 

g 



-l 

3 



l 

/ {G . } <G . >dn (9 _ } 
0 1 J J 



d-p)p c c c 



1 

/ {G . } <G . >dr| { 9 . } 
0 1 3 L 



{G i }L" 2 pk a <G i >'{9 2 } | J - L _2 pk a / {G.. } ' <G ^ > ' dr\ { 9 2 } (C.18) 

1 

- pp c uL / {G. }<G .> ' dri { 9 » } + hz /{G. } <G . >dr| { 9 } 

a a ~ 1 J Z 1 J -L 



-hz / {G . } <G . >dp { 9 - } + uL -1 (pP) ' / {G. }dp 

0 1 3 0 1 



= p p 



a 



c 

a 



1 

/ {G. }<G.>dn(0 9 } 
1 j ^ 
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(C . 19 ) 




The first term in each of the above expressions is a boundary 
term which permits incorporation of natural boundary condi- 
tions. Implementation of the boundary conditions is presented 
in Section C.l.b. The coefficients in equations C.17, C.18, 
and C.19 are comprised of variable dependent properties, and 
were taken as the average value of the properties over an 
element. In the limit, as the elements get smaller (i.e., 
n ->• oo ) , the average values of the coefficients converge to 
the exact values. 

Inspection of expressions C.17, C.18 and C.19 shows 
the four operators. 



/ { G ± } ' <Gj > 'dn 



(C . 20 ) 



/ {G i }<G j >'dn 



(C. 21) 



/ {G . } <G . >dr) 
; i 3 



(C . 22 ) 



/ {G^dn 



(C . 23) 
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To formulate these operators, the global shape function. 



G^, was defined on the local level by. 



G. = g 



(i-1) (Qv (i) 
l (+) 9? 



(C . 24 ) 



where g^ and g^ were defined by. 



( e ) 
g l 



(1 " l ] 
x e 



for E, in element (e) 



for E, not in element (e) 



(C . 25) 



JL 

JL 



for E, in element (e) 



.(e) 



(C . 26 ) 



for E not i n element (e) 



and is the length of the eth element. The (+) notation 
in expression C.24 means that G^ is the union of g|^ ^ and 
g^^ . The local shape functions (i.e., the elements) have 
the following properties. 



i 

0 

(i) / g| j) g/ m) =0 if j * m (C.27) 

0 1 K 

1 if i = j 

(ii) g i (e) (n j ) = 6 ^ = (C . 28 ) 

0 if i M 
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Having defined the local shape functions, the ele- 
mental matrix operators contributing to the global matrix 
operators C.20 through C.23 are. 



1 

/ {G . } ' <G . > ' dn 
0 1 3 





0 



1 

/ {G ± } <Gj > ' 




-1 

-1 



1 

1 



0 



/ 1 {G i }<G j 



> -> 




1 

2 



1 



/ (G } 

0 1 





(C. 29) 



(C • 30 ) 



(C . 31 ) 



(C. 32) 



The derivations of these operators are presented in Section 
C.l.d of this appendix. 

b. Implementation of Boundary Conditions 

Having formulated the system matrices for the field 
equations, treatment of the boundary conditions will now be 
discussed. Each field equation is considered individually. 

1 . Porous Solid Transfer Equation 

The third set of porous solid heat transfer 
boundary conditions (i.e., expressions III. 57 and III. 58) will 
only be considered since the first and second sets are sub- 
sets of the third. The third set of boundary conditions for 
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the porous 


solid 


are , 






L'^VkJ 


a* £ 

II 


~4 ~4 

h, (T -T ) + as (T -T ) 
1 c 00 c 00 


n = o 


(C.33) 




II 

& 


~4 ~4 

-h_ (T -T ) - as (T -T ) 
2 c 00 c °° 


13 

II 


(C.34) 



Since the first term in expression C.17 is. 



{G. }L -2 (1-p) (k +k ) <G . > 1 { 0 , } | (C . 35) 

i e r i i o 



or in analogous form, 

L~ 2 ( 1-p) (k e +k r ) 

Natural boundary conditions, C.33 and C.34, may be directly 
substituted into equation C.36. The response dependent 
parameters, h^, , and T c , changing with time, are evaluated 

continuously within the integration routine. Thus, the boun- 
dary conditions are incorporated in the system matrices as 
follows . 

(1) -L~^(l-p)h^: added to the stiffness matrix A(t) 

at location A^ ^ ~ 

(2) L l(l-p) [h.T^-as (T^-T^) ] : added to the excitation 

vector F(t) at location 



96 . 



9n 



(C.36) 
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(3) -L ■*‘(l-p)h_: added to the stiffness matrix A(t) 

at location A- ~ ~ 

3n-2 , 3n-2 



(4) L -1 (l-p) [h 2 T oo -oe (T^-T^) ] : 



added to the excitation 
vector F(t) at location 
F 

r 3n-2 



2. Air Heat Transfer Equation 

The first and second set of boundary conditions 
for the air heat transfer equation are. 



i 3T 
L _1 k —S. 
a 3n 



p c u(T -T ) 
M a a a 



n = 0 



(C . 37) 



3T 

a 

an 



= o 



n = l 



(C. 38) 



Since these are both natural boundary conditions, they are 
substututed for the first term of expression C.18 at n = 0 
and n = 1, respectively. The time dependent properties and 
parameters in the coefficient are evaluated continuously 
within the integration routine. The boundary conditions are 
implemented by adding. 



Cl) -p p c uL to the stiffness matrix A(t) at location 

a a A = 

A 2,2 



(2) p p c uL^T : to the excitation vector F(t) at 



a a 



location F. 



The third set of air heat transfer boundary conditions are. 



T = T 
a °° 



n = 0 



(C. 39) 
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0 



n 



1 



(C . 40) 



3T 



9n 



The essential boundary condition at n = 0 is imposed in the 
Galerkin equation as follows. The A~ . row of the A(t) 
matrix, the B 9 • row and the B. 9 column of the B(t) matrix, 
and the F ^ location of the excitation vector, F(t), are all 
set equal to zero. The H >2 2 l° cat i° n °f the B(t) matrix 
is then set equal to one. 

3. Oxygen Transport Equation 

For the oxygen diffusion equation, the first and 
second set of boundary conditions are, 

L_1 V e Hr = U(<|) " O n = ° (C. 41) 

= 0 n = 1 (C. 42) 



Since these are natural boundary conditions, they were sub- 
stituted for the first term in expression C.19. The proper- 
ties are evaluated continuously within the integration routine. 
The boundary conditions are implemented by adding. 



(1) -p u L 1 : to the stiffness matrix A(t) at location 

A 3,3 



(2) p u L ^ ip : to the excitation vector F(t) at location 
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The third set of boundary conditions for the oxygen diffusion 
equation are. 



<f> = n = 0 (C. 43) 

1^-0 n = 1 (C. 44 ) 

The essential boundary condition at n = 0 is imposed in the 
Galerkin equation as follows. The A 7 . row of the A(t) 
matrix, the B 7 . row and the B. column of the B(t) matrix, 
and the F, location of the excitation vector, F(t), are all 

-J — 

set equal to zero. The B^ ^ location of the B(t) matrix is 
then set equal to one. 

4. Surface Recession Problem 

As discussed previously, the oxygen diffusion 
equation is eliminated during the surface recession phase of 
a problem. By doing this, a reordering of the system matrices 
must take place, including the locations for applying the 
boundary conditions. The boundary conditions of the particle 
and air temperature equations are implemented the same as above 
except the indices identifying the matrix locations of the 
form, 3n-i, are changed to 2n-(i-l). 

c . Treatment of the Reaction Rate Term 

An exponential reaction rate term appears in both the 
porous solid heat transfer and the oxygen diffusion equation. 
The modified Gear integration routine requires calculating or 
approximating the Jacobian matrix from the system matrices. 
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A(t) and B(t), In order to include the reaction rate term 
in the Jacobian matrix, thus improving the efficiency of the 
integration routine, the exponential terms were treated as 
follows. The reaction rate terms were multiplied and divided 
by the oxygen concentration (only the heat generation term 
is shown, oxygen consumption is treated identically) , 



R 

g 




l 



i=l 



V 3i 



(C. 45) 



Applying the operation of expression C.14, and moving the 
reaction rate term divided by the oxygen concentration out- 
side the integral, the Galerkin operator for the heat 
generation becomes, 

1 1 

R g z9 3 / <G j >dri ( 9 3 > (C. 46) 

The coefficient of expression C.46 is treated as a time 
dependent parameter. The Galerkin operators are distributed 
into the stiffness matrix at locations, 3 n-2, 3 n for heat 
generation, and 3 n, 3 n for oxygen consumption. During the 
surface recession phase, reaction rate is treated in a 
different manner as discussed in Section III.G. 
d. Derivation of the FEM Operators 

In the section on the finite element formulation, 
the following four differential operators were identified. 
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(C. 47) 



1 

/ {G.'}<G!>dn 
0 3 


(C. 47) 


1 

/ {G . }<G ! >dn 
0 3 


(C. 48) 


1 

/ {G . }<G . >dri 
0 1 3 


(C.49) 


1 

/ {G. }dn 
0 1 


(C . 50 ) 



where the are the global basis functions. These operators 
are constructed on the element level by introducing the 
corresponding element basis functions, g^. The global and 
element basis functions are related by. 



. _ (i-1) ^ _ (i) 

G i - & ?2 

where g^ and g^ are defined by 


(C . 51) 


(1 - -f-) for £ in element (e) 

(e) _ S 
*1 

0 for £ not in element (e) 


(C . 52 ) 
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for 5 in element (e) 



_L 

( e) = 

g 2 ~ (C . 53 ) 

0 for C not in element (e) 



and l is the length of the (e)th element. 

The derivation of the local elemental matrices 
according to the Galerkin method for the global operations 
proceeds as follows; 

For operator C.47, 



global 



local 



1 

/ {G'}<G!>dn 
0 X 3 




(C . 54 ) 



Noting that , 




(C.55) 



94 




the elemental matrix becomes 



(C . 56 ) 



£ 

e 

/ 


1 

to 


-(-L) 2 - 
v £ ’ 
e 


<35 - T~ 


~ i 


-r 


0 


L-<f > 2 

e 


n> 

to 

1 


e 


_-i 


i_ 



(C . 57) 
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For operator C.48, 



global 



local 



/ {G i }<G!>dn 




<g{ g*>d£] (C . 58) 



Substituting in the local shape functions gives. 



(1 - j-) 

e 



<f> 

e 



< ( - <f)>d 5 

e e 



(C. 59) 



Carrying out the operations gives 



( 1 . 5 > , 1 € n 




'-1 1" 


e e e * e 


i 




O 0 


dC - 2 




r 2 £ 2 

L- ( fr) <£t> J 




_-l 1_ 



21 



(C . 60 ) 



For operator C.49, 
global 



local 



/ {G . } <G . >dn 

. i"i 



Ut / 

0 




<g x g 2 >dC] (C . 61) 



**2 
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Substituting in for the local shape functions 




(C. 62) 



the elemental matrix becomes 





1 - 
2 



(C. 63) 



1 J 



For operator C.50, 



global 



local 



/ 

0 



1 



{G i >dn 




(C . 64 ) 



Substituting in the local shape function, and integrating. 



the expression becomes 



l 



e 



/ 

0 





(C . 65) 
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This last operator is used for the excitation vector as 
described in the FEM formulation. 

2. SHOOTING METHOD 

The solution of equation III. 12 is obtained by the 
shooting method, a general discussion of the shooting method 
with examples is presented by Gerald [55] . The method is 
based on iterative solutions of equation III. 12 until a solu- 
tion is reached which also satisfies the boundary conditions, 
P(0) = and P(l) = P L . Along with the boundary condition 
at n =0, an initial estimate to dP/dri is specified, and 
equation III. 12 is integrated from ri = 0 to ri = 1 using 
Euler's method. An approximate Newton's method is then used 
as follows. 



dP 

dn 



i 



£| - p <!> 
an | j.! 



dP 

dn 



dP 

dn 



i-2 



i-1 P(l), - P(l) 



(C . 66 ) 



i-2 



to provide a better approximation of dP/dn . The procedure is 
repeated until the solution has converged. A solution satis- 
fying equation III. 12 and its associated boundary conditions 
is calculated at each time step. The current values of the 
properties are used and 9p /3t is approximated by linear 

cl 

t 

interpolation using the current and past values of p . The 

cl 

3p /dt term was neglected for the first two integrations of a 
problem. The reason for this is associated with the difficulty 
of specifying a reasonable set of initial conditions for the 
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temperature and oxygen diffusion equations. When equation 
III. 12 was treated as an initial condition problem, the 
pressure gradient in a small region of the porous medium would 
approach zero for some initial conditions of temperature and 
oxygen concentration (the starting heat generation rate for 
those initial conditions was large at time, t = 0) . The 
pressure gradient approaching zero resulted in numerical 
instabilities for the temperature and oxygen diffusion equations. 
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